diff --git a/include/base/dof_object.h b/include/base/dof_object.h
index 655dda52752..5f8669e2ca0 100644
--- a/include/base/dof_object.h
+++ b/include/base/dof_object.h
@@ -181,12 +181,14 @@ class DofObject : public ReferenceCountedObject<DofObject>
    */
   unique_id_type unique_id () const;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The globally \p unique_id for this \p DofObject as a
    * writable reference.  Deprecated; use the API taking an input
    * instead.
    */
   unique_id_type & set_unique_id ();
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Sets the \p unique_id for this \p DofObject
@@ -851,6 +853,7 @@ unique_id_type DofObject::unique_id () const
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 inline
 unique_id_type & DofObject::set_unique_id ()
 {
@@ -861,6 +864,7 @@ unique_id_type & DofObject::set_unique_id ()
   libmesh_not_implemented();
 #endif
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/include/fe/fe.h b/include/fe/fe.h
index 1b9873e6cbe..8922fd47011 100644
--- a/include/fe/fe.h
+++ b/include/fe/fe.h
@@ -414,6 +414,11 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
    * a finite element of type \p t and approximation order \p o.
    *
    * On a p-refined element, \p o should be the total order of the element.
+   *
+   * This method does not support all finite element types; e.g. for an
+   * arbitrary polygon or polyhedron type the number of shape
+   * functions may depend on an individual element and not just its
+   * type.
    */
   static unsigned int n_shape_functions (const ElemType t,
                                          const Order o)
@@ -424,29 +429,72 @@ class FE : public FEGenericBase<typename FEOutputType<T>::type>
    * finite element.
    *
    * On a p-refined element, \p o should be the total order of the element.
+   *
+   * This method does not support all finite element types; e.g. for an
+   * arbitrary polygon or polyhedron type the number of shape
+   * functions may depend on an individual element and not just its
+   * type.
    */
   static unsigned int n_dofs(const ElemType t,
                              const Order o);
 
+  /**
+   * \returns The number of shape functions associated with this
+   * finite element.
+   *
+   * On a p-refined element, \p o should be the total order of the element.
+   *
+   * \p e should only be a null pointer if using a FE family like
+   * SCALAR that has degrees of freedom independent of any element.
+   */
+  static unsigned int n_dofs(const Elem * e,
+                             const Order o);
+
   /**
    * \returns The number of dofs at node \p n for a finite element
    * of type \p t and order \p o.
    *
    * On a p-refined element, \p o should be the total order of the element.
+   *
+   * This method does not support all finite element types; e.g. for an
+   * arbitrary polygon or polyhedron type the meaning of a node index
+   * \p n may depend on an individual element and not just its type.
    */
   static unsigned int n_dofs_at_node(const ElemType t,
                                      const Order o,
                                      const unsigned int n);
 
+  /**
+   * \returns The number of dofs at node \p n for a finite element
+   * of type \p t and order \p o.
+   *
+   * On a p-refined element, \p o should be the total order of the element.
+   */
+  static unsigned int n_dofs_at_node(const Elem & e,
+                                     const Order o,
+                                     const unsigned int n);
+
   /**
    * \returns The number of dofs interior to the element,
    * not associated with any interior nodes.
    *
    * On a p-refined element, \p o should be the total order of the element.
+   *
+   * This method may not support all finite element types, e.g. higher
+   * order polygons or polyhedra may differ from element to element.
    */
   static unsigned int n_dofs_per_elem(const ElemType t,
                                       const Order o);
 
+  /**
+   * \returns The number of dofs interior to the element,
+   * not associated with any interior nodes.
+   *
+   * On a p-refined element, \p o should be the total order of the element.
+   */
+  static unsigned int n_dofs_per_elem(const Elem & e,
+                                      const Order o);
+
   /**
    * \returns The continuity level of the finite element.
    */
@@ -1528,6 +1576,8 @@ void lagrange_nodal_soln(const Elem * elem,
  */
 unsigned int monomial_n_dofs(const ElemType t, const Order o);
 
+unsigned int monomial_n_dofs(const Elem * e, const Order o);
+
 /**
  * Helper functions for rational basis functions.
  */
@@ -1582,6 +1632,73 @@ void rational_all_shape_derivs (const Elem & elem,
 
 } // namespace libMesh
 
+
+// Full specialization of all n_dofs type functions, for every
+// dimension, with both original ElemType and new Elem signatures
+#define LIBMESH_DEFAULT_NDOFS(MyType) \
+template <> unsigned int FE<0,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
+template <> unsigned int FE<1,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
+template <> unsigned int FE<2,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
+template <> unsigned int FE<3,MyType>::n_dofs(const ElemType t, const Order o) { return MyType##_n_dofs(t, o); } \
+ \
+template <> unsigned int FE<0,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
+template <> unsigned int FE<1,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
+template <> unsigned int FE<2,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
+template <> unsigned int FE<3,MyType>::n_dofs(const Elem * e, const Order o) { return MyType##_n_dofs(e, o); } \
+ \
+template <> unsigned int FE<0,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
+template <> unsigned int FE<1,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
+template <> unsigned int FE<2,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
+template <> unsigned int FE<3,MyType>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(t, o, n); } \
+ \
+template <> unsigned int FE<0,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
+template <> unsigned int FE<1,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
+template <> unsigned int FE<2,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
+template <> unsigned int FE<3,MyType>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return MyType##_n_dofs_at_node(e, o, n); } \
+ \
+template <> unsigned int FE<0,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
+template <> unsigned int FE<1,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
+template <> unsigned int FE<2,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
+template <> unsigned int FE<3,MyType>::n_dofs_per_elem(const ElemType t, const Order o) { return MyType##_n_dofs_per_elem(t, o); } \
+ \
+template <> unsigned int FE<0,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
+template <> unsigned int FE<1,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
+template <> unsigned int FE<2,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); } \
+template <> unsigned int FE<3,MyType>::n_dofs_per_elem(const Elem & e, const Order o) { return MyType##_n_dofs_per_elem(e, o); }
+
+
+#define LIBMESH_DEFAULT_VEC_NDOFS(MyType) \
+template <> unsigned int FE<0,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,MyType>::n_dofs(t, o); } \
+template <> unsigned int FE<1,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,MyType>::n_dofs(t, o); } \
+template <> unsigned int FE<2,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,MyType>::n_dofs(t, o); } \
+template <> unsigned int FE<3,MyType##_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,MyType>::n_dofs(t, o); } \
+ \
+template <> unsigned int FE<0,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,MyType>::n_dofs(e, o); } \
+template <> unsigned int FE<1,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,MyType>::n_dofs(e, o); } \
+template <> unsigned int FE<2,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,MyType>::n_dofs(e, o); } \
+template <> unsigned int FE<3,MyType##_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,MyType>::n_dofs(e, o); } \
+ \
+template <> unsigned int FE<0,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,MyType>::n_dofs_at_node(t, o, n); } \
+template <> unsigned int FE<1,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,MyType>::n_dofs_at_node(t, o, n); } \
+template <> unsigned int FE<2,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 2*FE<2,MyType>::n_dofs_at_node(t, o, n); } \
+template <> unsigned int FE<3,MyType##_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<3,MyType>::n_dofs_at_node(t, o, n); } \
+ \
+template <> unsigned int FE<0,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,MyType>::n_dofs_at_node(e.type(), o, n); } \
+template <> unsigned int FE<1,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,MyType>::n_dofs_at_node(e.type(), o, n); } \
+template <> unsigned int FE<2,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 2*FE<2,MyType>::n_dofs_at_node(e.type(), o, n); } \
+template <> unsigned int FE<3,MyType##_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 3*FE<3,MyType>::n_dofs_at_node(e.type(), o, n); } \
+ \
+template <> unsigned int FE<0,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,MyType>::n_dofs_per_elem(t, o); } \
+template <> unsigned int FE<1,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,MyType>::n_dofs_per_elem(t, o); } \
+template <> unsigned int FE<2,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 2*FE<2,MyType>::n_dofs_per_elem(t, o); } \
+template <> unsigned int FE<3,MyType##_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 3*FE<3,MyType>::n_dofs_per_elem(t, o); } \
+ \
+template <> unsigned int FE<0,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,MyType>::n_dofs_per_elem(e.type(), o); } \
+template <> unsigned int FE<1,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,MyType>::n_dofs_per_elem(e.type(), o); } \
+template <> unsigned int FE<2,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 2*FE<2,MyType>::n_dofs_per_elem(e.type(), o); } \
+template <> unsigned int FE<3,MyType##_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 3*FE<3,MyType>::n_dofs_per_elem(e.type(), o); }
+
+
 #define LIBMESH_DEFAULT_VECTORIZED_FE(MyDim, MyType) \
 template<>                                           \
 void FE<MyDim,MyType>::all_shapes                    \
diff --git a/include/fe/fe_abstract.h b/include/fe/fe_abstract.h
index 984e56e10fc..ad8d8b7d53d 100644
--- a/include/fe/fe_abstract.h
+++ b/include/fe/fe_abstract.h
@@ -201,16 +201,24 @@ class FEAbstract : public ReferenceCountedObject<FEAbstract>
                          const std::vector<Point> & reference_side_points,
                          std::vector<Point> &       reference_points) = 0;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns \p true if the point p is located on the reference element
    * for element type t, false otherwise.  Since we are doing floating
    * point comparisons here the parameter \p eps can be specified to
    * indicate a tolerance.  For example, \f$ x \le 1 \f$  becomes
    * \f$ x \le 1 + \epsilon \f$.
+   *
+   * \deprecated This method overload does not support all finite
+   * element types; e.g. the reference element for an arbitrary
+   * polygon or polyhedron type may differ from element to element.
+   * Use \p Elem::on_reference_element() instead.
    */
   static bool on_reference_element(const Point & p,
                                    const ElemType t,
                                    const Real eps = TOLERANCE);
+#endif // LIBMESH_ENABLE_DEPRECATED
+
   /**
    * \returns The reference space coordinates of \p nodes based on the
    * element type.
diff --git a/include/fe/fe_interface.h b/include/fe/fe_interface.h
index a747ce6401c..622f2d7f5ef 100644
--- a/include/fe/fe_interface.h
+++ b/include/fe/fe_interface.h
@@ -78,21 +78,21 @@ class FEInterface
    */
   virtual ~FEInterface() = default;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
-   * \returns The number of shape functions associated with this
-   * finite element of type \p fe_t.
-   * Automatically decides which finite element class to use.
-   *
-   * On a p-refined element, \p fe_t.order should be the total order of the element.
-   *
    * \deprecated Call the version of this function taking an Elem* instead.
    */
   static unsigned int n_shape_functions(const unsigned int dim,
                                         const FEType & fe_t,
                                         const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
-   * Non-deprecated version of function above.
+   * \returns The number of shape functions associated with this
+   * finite element \p elem of type \p fe_t.
+   * Automatically decides which finite element class to use.
+   *
+   * On a p-refined element, \p fe_t.order should be the total order of the element.
    */
   static unsigned int n_shape_functions(const FEType & fe_t,
                                         const Elem * elem,
@@ -106,13 +106,8 @@ class FEInterface
                                         const int extra_order,
                                         const Elem * elem);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
-   * \returns The number of shape functions associated with this
-   * finite element.
-   * Automatically decides which finite element class to use.
-   *
-   * On a p-refined element, \p fe_t.order should be the total order of the element.
-   *
    * \deprecated Use n_dofs(const FEType &, Elem*) or n_dofs(const FEType &, int, Elem*) instead.
    */
   static unsigned int n_dofs(const unsigned int dim,
@@ -131,6 +126,7 @@ class FEInterface
   static unsigned int n_dofs(const unsigned int dim,
                              const FEType & fe_t,
                              const Elem * elem);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * \returns The number of DOFs for \p elem for finite element type \p fe_t
@@ -152,6 +148,11 @@ class FEInterface
                              int extra_order,
                              const Elem * elem);
 
+  typedef unsigned int (*n_dofs_at_node_ptr) (const ElemType,
+                                              const Order,
+                                              const unsigned int);
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The number of dofs at node n for a finite element
    * of type \p fe_t.
@@ -168,14 +169,7 @@ class FEInterface
                                      const ElemType t,
                                      const unsigned int n);
 
-  typedef unsigned int (*n_dofs_at_node_ptr) (const ElemType,
-                                              const Order,
-                                              const unsigned int);
-
   /**
-   * \returns A function which evaluates n_dofs_at_node for the
-   * requested FE type and dimension.
-   *
    * \deprecated Use the version of this function that takes an Elem*
    * for consistency. The behavior is otherwise exactly the same,
    * since this function does not depend on the Elem::p_level().
@@ -183,9 +177,11 @@ class FEInterface
   static n_dofs_at_node_ptr
   n_dofs_at_node_function(const unsigned int dim,
                           const FEType & fe_t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
-   * Non-deprecated version of function above.
+   * \returns A function which evaluates n_dofs_at_node for the
+   * requested FE type and element.
    */
   static n_dofs_at_node_ptr
   n_dofs_at_node_function(const FEType & fe_t,
@@ -211,21 +207,21 @@ class FEInterface
                                      const Elem * elem,
                                      const unsigned int n);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
-   * \returns The number of dofs interior to the element,
-   * not associated with any interior nodes.
-   * Automatically decides which finite element class to use.
-   *
-   * On a p-refined element, \p fe_t.order should be the total order of the element.
-   *
    * \deprecated Call the version of this function that takes an Elem* instead.
    */
   static unsigned int n_dofs_per_elem(const unsigned int dim,
                                       const FEType & fe_t,
                                       const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
-   * The non-deprecated version of the function above.
+   * \returns The number of dofs interior to the element,
+   * not associated with any interior nodes.
+   * Automatically decides which finite element class to use.
+   *
+   * On a p-refined element, \p fe_t.order should be the total order of the element.
    */
   static unsigned int n_dofs_per_elem(const FEType & fe_t,
                                       const Elem * elem,
@@ -311,6 +307,7 @@ class FEInterface
                               std::vector<Number> & nodal_soln,
                               const bool add_p_level = true);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * This is now deprecated; use FEMap::map instead.
    */
@@ -347,10 +344,16 @@ class FEInterface
    * Since we are doing floating point comparisons, the parameter
    * \p eps can be specified to indicate a tolerance.  For example,
    * \f$ \xi \le 1 \f$ becomes \f$ \xi \le 1 + \epsilon \f$.
+   *
+   * \deprecated This method overload does not support all finite
+   * element types; e.g. the reference element for an arbitrary
+   * polygon or polyhedron type may differ from element to element.
+   * Use \p Elem::on_reference_element() instead.
    */
   static bool on_reference_element(const Point & p,
                                    const ElemType t,
                                    const Real eps=TOLERANCE);
+
   /**
    * \returns The value of the \f$ i^{th} \f$ shape function at
    * point \p p. This method allows you to specify the dimension,
@@ -388,8 +391,12 @@ class FEInterface
                     const Elem * elem,
                     const unsigned int i,
                     const Point & p);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
+   * \returns The value of the \f$ i^{th} \f$ shape function at
+   * point \p p.
+   *
    * Non-deprecated version of the shape() function. The
    * Elem::p_level() is accounted for internally if \p add_p_level
    */
@@ -400,6 +407,9 @@ class FEInterface
                     const bool add_p_level = true);
 
   /**
+   * \returns The value of the \f$ i^{th} \f$ shape function at
+   * point \p p.
+   *
    * Non-deprecated version of the shape() function. The
    * Elem::p_level() is ignored and extra_order is used instead.
    */
@@ -409,6 +419,7 @@ class FEInterface
                     const unsigned int i,
                     const Point & p);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The value of the \f$ i^{th} \f$ shape function at
    * point \p p. This method allows you to specify the dimension,
@@ -450,8 +461,14 @@ class FEInterface
                     const unsigned int i,
                     const Point & p,
                     OutputType & phi);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
+   * \returns The value of the \f$ i^{th} \f$ shape function at
+   * point \p p. This method allows you to specify the dimension,
+   * element type, and order directly. Automatically passes the
+   * request to the appropriate *scalar* finite element class member.
+   *
    * Non-deprecated version of templated shape() function that
    * accounts for Elem::p_level() internally.
    */
@@ -463,6 +480,11 @@ class FEInterface
                     OutputType & phi);
 
   /**
+   * \returns The value of the \f$ i^{th} \f$ shape function at
+   * point \p p. This method allows you to specify the dimension,
+   * element type, and order directly. Automatically passes the
+   * request to the appropriate *scalar* finite element class member.
+   *
    * Non-deprecated version of templated shape() function that ignores
    * Elem::p_level() and instead uses extra_order internally.
    */
@@ -520,24 +542,25 @@ class FEInterface
                              const Point & p,
                              const bool add_p_level);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
-   * \returns A function which evaluates shape for the
-   * requested FE type and dimension.
-   *
    * \deprecated Call the version of this function taking an Elem* instead.
    */
   static shape_ptr
   shape_function(const unsigned int dim,
                  const FEType & fe_t,
                  const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
-   * Non-deprecated version of function above.
+   * \returns A function which evaluates shape for the
+   * requested FE type and element.
    */
   static shape_ptr
   shape_function(const FEType & fe_t,
                  const Elem * elem);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The \f$ j^{th} \f$ coordinate of the gradient of
    * the \f$ i^{th} \f$ shape function at point \p p.
@@ -558,15 +581,6 @@ class FEInterface
                           const Point & p);
 
   /**
-   * \returns The \f$ j^{th} \f$ coordinate of the gradient of
-   * the \f$ i^{th} \f$ shape function at point \p p.
-   * This method allows you to specify the dimension,
-   * element type, and order directly. Automatically passes the
-   * request to the appropriate *scalar* finite element class member.
-   *
-   * \note On a p-refined element, \p fe_t.order should be the total
-   * order of the element.
-   *
    * \deprecated Call the version of this function taking an Elem* instead.
    */
   static Real shape_deriv (const unsigned int dim,
@@ -575,9 +589,17 @@ class FEInterface
                            const unsigned int i,
                            const unsigned int j,
                            const Point & p);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
-   * Non-deprecated version of function above.
+   * \returns The \f$ j^{th} \f$ coordinate of the gradient of
+   * the \f$ i^{th} \f$ shape function at point \p p.
+   * This method allows you to specify the dimension,
+   * element, and order directly. Automatically passes the
+   * request to the appropriate *scalar* finite element class member.
+   *
+   * \note On a p-refined element, \p fe_t.order should be the total
+   * order of the element.
    */
   static Real shape_deriv(const FEType & fe_t,
                           const Elem * elem,
@@ -647,6 +669,8 @@ class FEInterface
                        const Elem * elem);
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The second \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
    * shape function at the point \p p.
@@ -676,6 +700,18 @@ class FEInterface
                                  const unsigned int j,
                                  const Point & p);
 
+  /**
+   * \deprecated Call version of this function which does not require
+   * \p dim and takes an Elem * instead.
+   */
+  static Real shape_second_deriv (const unsigned int dim,
+                                  const FEType & fe_t,
+                                  const Elem *elem,
+                                  const unsigned int i,
+                                  const unsigned int j,
+                                  const Point & p);
+#endif // LIBMESH_ENABLE_DEPRECATED
+
   /**
    * \returns The second \f$ j^{th} \f$ derivative of the \f$ i^{th} \f$
    * shape function at the point \p p.
@@ -690,19 +726,6 @@ class FEInterface
    *
    * \note On a p-refined element, \p fe_t.order should be the total
    * order of the element.
-   *
-   * \deprecated Call version of this function which does not require
-   * \p dim and takes an Elem * instead.
-   */
-  static Real shape_second_deriv (const unsigned int dim,
-                                  const FEType & fe_t,
-                                  const Elem *elem,
-                                  const unsigned int i,
-                                  const unsigned int j,
-                                  const Point & p);
-
-  /**
-   * Non-deprecated version of function above.
    */
   static Real shape_second_deriv(const FEType & fe_t,
                                  const Elem * elem,
@@ -732,19 +755,19 @@ class FEInterface
                                           const Point & p,
                                           const bool add_p_level);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
-   * \returns A function which evaluates shape for the
-   * requested FE type and dimension.
-   *
    * \deprecated Call the version of this function that takes an Elem * instead.
    */
   static shape_second_deriv_ptr
   shape_second_deriv_function(const unsigned int dim,
                               const FEType & fe_t,
                               const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
-   * Non-deprecated version of the function above.
+   * \returns A function which evaluates shape for the
+   * requested FE type and dimension.
    */
   static shape_second_deriv_ptr
   shape_second_deriv_function(const FEType & fe_t,
@@ -862,6 +885,7 @@ class FEInterface
    * the calls to \p FE and \p InfFE.
    */
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \deprecated Call the version of ifem_n_shape_functions() which
    * takes a pointer-to-Elem instead.
@@ -869,10 +893,12 @@ class FEInterface
   static unsigned int ifem_n_shape_functions(const unsigned int dim,
                                              const FEType & fe_t,
                                              const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   static unsigned int ifem_n_shape_functions(const FEType & fe_t,
                                              const Elem * elem);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \deprecated Call the version of ifem_n_dofs() which takes a
    * pointer-to-Elem instead.
@@ -880,10 +906,12 @@ class FEInterface
   static unsigned int ifem_n_dofs(const unsigned int dim,
                                   const FEType & fe_t,
                                   const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   static unsigned int ifem_n_dofs(const FEType & fe_t,
                                   const Elem * elem);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \deprecated Call the version of ifem_n_dofs_at_node() which takes
    * a pointer-to-Elem instead.
@@ -892,11 +920,13 @@ class FEInterface
                                           const FEType & fe_t,
                                           const ElemType t,
                                           const unsigned int n);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   static unsigned int ifem_n_dofs_at_node(const FEType & fe_t,
                                           const Elem * elem,
                                           const unsigned int n);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \deprecated Call the version of ifem_n_dofs_per_elem() which
    * takes a pointer-to-Elem instead.
@@ -904,6 +934,7 @@ class FEInterface
   static unsigned int ifem_n_dofs_per_elem(const unsigned int dim,
                                            const FEType & fe_t,
                                            const ElemType t);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   static unsigned int ifem_n_dofs_per_elem(const FEType & fe_t,
                                            const Elem * elem);
@@ -934,7 +965,11 @@ class FEInterface
                                 const Real tolerance = TOLERANCE,
                                 const bool secure = true);
 
-
+#ifdef LIBMESH_ENABLE_DEPRECATED
+  /*
+   * \deprecated This method overload may not support all finite
+   * element types.  Use \p Elem::on_reference_element() instead.
+   */
   static bool ifem_on_reference_element(const Point & p,
                                         const ElemType t,
                                         const Real eps);
@@ -958,12 +993,14 @@ class FEInterface
                          const Elem * elem,
                          const unsigned int i,
                          const Point & p);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   static Real ifem_shape(const FEType & fe_t,
                          const Elem * t,
                          const unsigned int i,
                          const Point & p);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \deprecated Call version that takes a pointer-to-Elem and does
    * not require an explicit dim parameter instead.
@@ -985,6 +1022,7 @@ class FEInterface
                                const unsigned int i,
                                const unsigned int j,
                                const Point & p);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   static Real ifem_shape_deriv(const FEType & fe_t,
                                const Elem * elem,
diff --git a/include/fe/inf_fe.h b/include/fe/inf_fe.h
index e15f3b93b4d..46d555705a8 100644
--- a/include/fe/inf_fe.h
+++ b/include/fe/inf_fe.h
@@ -377,7 +377,13 @@ class InfFE : public FEBase
   /**
    * \returns The number of shape functions associated with
    * a finite element of type \p t and approximation order \p o.
-   *
+   */
+  static unsigned int n_shape_functions (const FEType & fet,
+                                         const Elem * inf_elem)
+  { return n_dofs(fet, inf_elem); }
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
+  /*
    * \deprecated Call the version of this function that takes a
    * pointer-to-Elem instead.
    */
@@ -385,40 +391,42 @@ class InfFE : public FEBase
                                          const ElemType t)
   { return n_dofs(fet, t); }
 
-  static unsigned int n_shape_functions (const FEType & fet,
-                                         const Elem * inf_elem)
-  { return n_dofs(fet, inf_elem); }
-
   /**
-   * \returns The number of shape functions associated with this
-   * infinite element.  Currently, we have \p o_radial+1 modes in
-   * radial direction, and \code FE<Dim-1,T>::n_dofs(...) \endcode
-   * in the base.
-   *
    * \deprecated Call the version of this function that takes an Elem*
    * instead for consistency with other FEInterface::n_dofs() methods.
    */
   static unsigned int n_dofs(const FEType & fet,
                              const ElemType inf_elem_type);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
+  /**
+   * \returns The number of shape functions associated with this
+   * infinite element.  Currently, we have \p o_radial+1 modes in
+   * radial direction, and \code FE<Dim-1,T>::n_dofs(...) \endcode
+   * in the base.
+   */
   static unsigned int n_dofs(const FEType & fet,
                              const Elem * inf_elem);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The number of dofs at infinite element node \p n
    * (not dof!) for an element of type \p t and order \p o.
-   *
-   * \deprecated Call the version of this function that takes an Elem*
-   * instead for consistency with other FEInterface::n_dofs() methods.
    */
   static unsigned int n_dofs_at_node(const FEType & fet,
                                      const ElemType inf_elem_type,
                                      const unsigned int n);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
+  /**
+   * \returns The number of dofs at infinite element node \p n
+   * (not dof!) for an element of type \p t and order \p o.
+   */
   static unsigned int n_dofs_at_node(const FEType & fet,
                                      const Elem * inf_elem,
                                      const unsigned int n);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \returns The number of dofs interior to the element,
    * not associated with any interior nodes.
@@ -428,7 +436,12 @@ class InfFE : public FEBase
    */
   static unsigned int n_dofs_per_elem(const FEType & fet,
                                       const ElemType inf_elem_type);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
+  /**
+   * \returns The number of dofs interior to the element,
+   * not associated with any interior nodes.
+   */
   static unsigned int n_dofs_per_elem(const FEType & fet,
                                       const Elem * inf_elem);
 
@@ -1020,12 +1033,8 @@ class InfFE : public FEBase
                                          unsigned int & base_node,
                                          unsigned int & radial_node);
 
-  /**
-   * Computes the indices of shape functions in the base \p base_shape and
-   * in radial direction \p radial_shape (0 in the base, \f$ \ge 1 \f$ further
-   * out) associated to the shape with global index \p i of an infinite element
-   * of type \p inf_elem_type.
-   *
+#ifdef LIBMESH_ENABLE_DEPRECATED
+  /*
    * \deprecated Call the version of this function that takes an Elem * instead.
    */
   static void compute_shape_indices (const FEType & fet,
@@ -1033,7 +1042,14 @@ class InfFE : public FEBase
                                      const unsigned int i,
                                      unsigned int & base_shape,
                                      unsigned int & radial_shape);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
+  /**
+   * Computes the indices of shape functions in the base \p base_shape and
+   * in radial direction \p radial_shape (0 in the base, \f$ \ge 1 \f$ further
+   * out) associated to the shape with global index \p i of an infinite element
+   * \p inf_elem.
+   */
   static void compute_shape_indices (const FEType & fet,
                                      const Elem * inf_elem,
                                      const unsigned int i,
diff --git a/include/geom/cell_hex.h b/include/geom/cell_hex.h
index da7ccc7d9a3..be54d95a9b5 100644
--- a/include/geom/cell_hex.h
+++ b/include/geom/cell_hex.h
@@ -195,6 +195,9 @@ class Hex : public Cell
    */
   static const unsigned int edge_sides_map[12][2];
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/cell_inf_hex.h b/include/geom/cell_inf_hex.h
index 78b1c3be6d0..8c66175ae35 100644
--- a/include/geom/cell_inf_hex.h
+++ b/include/geom/cell_inf_hex.h
@@ -226,6 +226,9 @@ class InfHex : public InfCell
    */
   static const unsigned int edge_sides_map[8][2];
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/cell_inf_prism.h b/include/geom/cell_inf_prism.h
index 444c7d71d3c..4b7bb870ad5 100644
--- a/include/geom/cell_inf_prism.h
+++ b/include/geom/cell_inf_prism.h
@@ -209,6 +209,9 @@ class InfPrism : public InfCell
    */
   static const unsigned int edge_sides_map[6][2];
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/cell_prism.h b/include/geom/cell_prism.h
index 992f8eabf4f..d27e7cad9ff 100644
--- a/include/geom/cell_prism.h
+++ b/include/geom/cell_prism.h
@@ -174,6 +174,9 @@ class Prism : public Cell
    */
   static const unsigned int edge_sides_map[9][2];
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/cell_pyramid.h b/include/geom/cell_pyramid.h
index e211e324d3d..e30c6f91609 100644
--- a/include/geom/cell_pyramid.h
+++ b/include/geom/cell_pyramid.h
@@ -198,6 +198,9 @@ class Pyramid : public Cell
    */
   static const unsigned int edge_sides_map[8][2];
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/cell_tet.h b/include/geom/cell_tet.h
index 88cf57b15a6..c8d8813d6c9 100644
--- a/include/geom/cell_tet.h
+++ b/include/geom/cell_tet.h
@@ -219,6 +219,9 @@ class Tet : public Cell
    */
   static const unsigned int edge_sides_map[6][2];
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/edge.h b/include/geom/edge.h
index 9531a6b436c..45a594c306c 100644
--- a/include/geom/edge.h
+++ b/include/geom/edge.h
@@ -234,6 +234,9 @@ class Edge : public Elem
 
   ElemType side_type (const unsigned int s) const override final;
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/elem.h b/include/geom/elem.h
index 94bd415387e..ed9a4473549 100644
--- a/include/geom/elem.h
+++ b/include/geom/elem.h
@@ -442,11 +442,13 @@ class Elem : public ReferenceCountedObject<Elem>,
   virtual unsigned int local_edge_node(unsigned int edge,
                                        unsigned int edge_node) const = 0;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * This function is deprecated, call local_side_node(side, side_node) instead.
    */
   unsigned int which_node_am_i(unsigned int side,
                                unsigned int side_node) const;
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * \returns \p true if a vertex of \p e is contained
@@ -945,14 +947,18 @@ class Elem : public ReferenceCountedObject<Elem>,
    */
   virtual Order default_side_order () const { return default_order(); }
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
-   * Calls Elem::vertex_average() for backwards compatibility. This
-   * method has now been deprecated, so it will be removed at some
-   * point in the future.  Calls to Elem::centroid() in user code
-   * should be updated to either call Elem::vertex_average() or
-   * Elem::true_centroid() on a case by case basis.
+   * Calls Elem::vertex_average() for backwards compatibility.
+   *
+   * \deprecated This method has now been deprecated, so it will be
+   * removed at some point in the future.  Calls to Elem::centroid()
+   * in user code should be updated to either call
+   * Elem::vertex_average() or Elem::true_centroid() on a case by case
+   * basis.
    */
   virtual Point centroid () const;
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * \returns The "true" geometric centroid of the element, c=(cx, cy,
@@ -1067,6 +1073,17 @@ class Elem : public ReferenceCountedObject<Elem>,
    */
   virtual bool contains_point (const Point & p, Real tol=TOLERANCE) const;
 
+  /**
+   * \returns \p true if the master-space point p is contained in the
+   * reference element corresponding to this element, false otherwise.
+   *
+   * Since we are doing floating point comparisons here the parameter
+   * \p eps can be specified to indicate a tolerance.  For example,
+   * \f$ x \le 1 \f$  becomes \f$ x \le 1 + \epsilon \f$.
+   */
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const = 0;
+
   /**
    * \returns \p true if this element is "close" to the point p, where
    * "close" is determined by the tolerance tol.
diff --git a/include/geom/face_inf_quad.h b/include/geom/face_inf_quad.h
index 9f063297cf2..a583e97133b 100644
--- a/include/geom/face_inf_quad.h
+++ b/include/geom/face_inf_quad.h
@@ -263,6 +263,9 @@ class InfQuad : public Elem
 
   virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int n) const override;
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/face_quad.h b/include/geom/face_quad.h
index 497d3197212..a3771fb3bb3 100644
--- a/include/geom/face_quad.h
+++ b/include/geom/face_quad.h
@@ -205,6 +205,9 @@ class Quad : public Face
 
   virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int n) const override;
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/face_tri.h b/include/geom/face_tri.h
index 74d4f0870d2..8fd722d7d5e 100644
--- a/include/geom/face_tri.h
+++ b/include/geom/face_tri.h
@@ -192,6 +192,9 @@ class Tri : public Face
 
   virtual std::vector<unsigned int> edges_adjacent_to_node(const unsigned int n) const override;
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/node_elem.h b/include/geom/node_elem.h
index ac622342a94..974ec8b61c1 100644
--- a/include/geom/node_elem.h
+++ b/include/geom/node_elem.h
@@ -315,6 +315,9 @@ class NodeElem : public Elem
    */
   virtual bool close_to_point(const Point & p, Real tol) const override;
 
+  virtual bool on_reference_element(const Point & p,
+                                    const Real eps = TOLERANCE) const override final;
+
 protected:
 
   /**
diff --git a/include/geom/remote_elem.h b/include/geom/remote_elem.h
index 6a98113a691..0cd0dc209f3 100644
--- a/include/geom/remote_elem.h
+++ b/include/geom/remote_elem.h
@@ -256,6 +256,10 @@ class RemoteElem : public Elem,
   virtual ElemType side_type (const unsigned int) const override
   { remote_elem_error("side_type"); return INVALID_ELEM; }
 
+  virtual bool on_reference_element(const Point & /*p*/,
+                                    const Real /*eps*/ = TOLERANCE) const override
+  { remote_elem_error("on_reference_element"); return false; }
+
 protected:
 
   /**
diff --git a/include/mesh/dyna_io.h b/include/mesh/dyna_io.h
index 3d0141487fa..6c58d93531f 100644
--- a/include/mesh/dyna_io.h
+++ b/include/mesh/dyna_io.h
@@ -99,12 +99,16 @@ class DynaIO : public MeshInput<MeshBase>
                               unsigned int sys_num,
                               unsigned int var_num);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * Removes any spline nodes (both NodeElem and Node), leaving only
    * the FE mesh generated from those splines.  Also removes node
    * constraints to the now-missing nodes.
+   *
+   * \p deprecated - use MeshTools::clear_spline_nodes(mesh) instead.
    */
   void clear_spline_nodes();
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * The integer type DYNA uses
diff --git a/include/mesh/mesh_base.h b/include/mesh/mesh_base.h
index b5769196611..a9f66363678 100644
--- a/include/mesh/mesh_base.h
+++ b/include/mesh/mesh_base.h
@@ -1143,8 +1143,10 @@ class MeshBase : public ParallelObject
    * will be deleted here - to keep those elements replicated during
    * preparation, set allow_remote_element_removal(false).
    */
+#ifdef LIBMESH_ENABLE_DEPRECATED
   void prepare_for_use (const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors);
   void prepare_for_use (const bool skip_renumber_nodes_and_elements);
+#endif // LIBMESH_ENABLE_DEPRECATED
   void prepare_for_use ();
 
   /**
@@ -1753,6 +1755,7 @@ class MeshBase : public ParallelObject
    */
   std::string get_local_constraints(bool print_nonlocal=false) const;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * \deprecated This method has ben replaced by \p cache_elem_data which
    * caches data in addition to elem dimensions (e.g. elem subdomain ids)
@@ -1761,6 +1764,7 @@ class MeshBase : public ParallelObject
    * be done manually by other classes after major mesh modifications.
    */
   void cache_elem_dims();
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /*
    * Search the mesh and cache data for the elements
diff --git a/include/mesh/mesh_communication.h b/include/mesh/mesh_communication.h
index 89d400809ca..c05b7b71f53 100644
--- a/include/mesh/mesh_communication.h
+++ b/include/mesh/mesh_communication.h
@@ -282,11 +282,13 @@ void connect_children(const MeshBase & mesh,
 // constraint rows, insert elements with the constraining nodes for
 // any constrained nodes in our set.
 //
-// This method is now deprecated, because it does not handle recursive
-// dependencies.  Use the new connect_element_dependencies method
-// instead.
+// \deprecated This method is now deprecated, because it does not
+// handle recursive dependencies.  Use the new
+// connect_element_dependencies method instead.
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void connect_families(connected_elem_set_type & connected_elements,
                       const MeshBase * mesh = nullptr);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 // Take a set of elements and create a set of connected nodes.
 void reconnect_nodes (connected_elem_set_type & connected_elements,
diff --git a/include/mesh/mesh_function.h b/include/mesh/mesh_function.h
index 5871146dc20..3e1ac6efe96 100644
--- a/include/mesh/mesh_function.h
+++ b/include/mesh/mesh_function.h
@@ -109,6 +109,7 @@ class MeshFunction : public FunctionBase<Number>,
    */
   virtual void init () override;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * The actual initialization process.  Takes an optional argument which
    * specifies the method to use when building a \p PointLocator
@@ -118,6 +119,7 @@ class MeshFunction : public FunctionBase<Number>,
    * instead call the version of init() taking no args.
    */
   void init (const Trees::BuildType point_locator_build_type);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Clears the function.
diff --git a/include/reduced_basis/rb_evaluation.h b/include/reduced_basis/rb_evaluation.h
index 136530a0dcc..bda15bf0c65 100644
--- a/include/reduced_basis/rb_evaluation.h
+++ b/include/reduced_basis/rb_evaluation.h
@@ -150,6 +150,7 @@ class RBEvaluation : public RBParametrized,
    */
   virtual Real residual_scaling_denom(Real alpha_LB);
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * Evaluate the dual norm of output \p n for the current parameters.
    *
@@ -161,6 +162,7 @@ class RBEvaluation : public RBParametrized,
    * theta values.
    */
   Real eval_output_dual_norm(unsigned int n, const RBParameters & mu);
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Evaluate the dual norm of output \p n for the current parameters,
diff --git a/include/reduced_basis/rb_parametrized.h b/include/reduced_basis/rb_parametrized.h
index 06aff38acaa..6e4c91660dc 100644
--- a/include/reduced_basis/rb_parametrized.h
+++ b/include/reduced_basis/rb_parametrized.h
@@ -99,6 +99,7 @@ class RBParametrized : public ReferenceCountedObject<RBParametrized>
    */
   unsigned int get_n_discrete_params() const;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * Get a set that stores the parameter names.
    *
@@ -108,6 +109,7 @@ class RBParametrized : public ReferenceCountedObject<RBParametrized>
    * loop over the keys directly.
    */
   std::set<std::string> get_parameter_names() const;
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Get the current parameters.
diff --git a/include/systems/diff_system.h b/include/systems/diff_system.h
index bb0e7a5602f..3de5fa627f3 100644
--- a/include/systems/diff_system.h
+++ b/include/systems/diff_system.h
@@ -201,11 +201,13 @@ class DifferentiableSystem : public ImplicitSystem,
   { this->_diff_physics.push(physics_in->clone_physics());
     this->_diff_physics.top()->init_physics(*this);}
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * Swap current physics object with external object.  This is
-   * deprecated.
+   * \deprecated for being too unsafe for easy use.
    */
   void swap_physics ( DifferentiablePhysics * & swap_physics );
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Push a clone of a new physics object onto our stack, overriding
diff --git a/include/systems/equation_systems.h b/include/systems/equation_systems.h
index 5900ef47529..cc7e90744b8 100644
--- a/include/systems/equation_systems.h
+++ b/include/systems/equation_systems.h
@@ -326,6 +326,7 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
   void get_vars_active_subdomains(const std::vector<std::string> & names,
                                   std::vector<std::set<subdomain_id_type>> & vars_active_subdomains) const;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * Retrieve the solution data for CONSTANT MONOMIALs and/or components of
    * CONSTANT MONOMIAL_VECs. If 'names' is populated, only the variables
@@ -337,6 +338,7 @@ class EquationSystems : public ReferenceCountedObject<EquationSystems>,
    */
   void get_solution (std::vector<Number> & soln,
                      std::vector<std::string> & names) const;
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Retrieve the solution data for CONSTANT MONOMIALs and/or components of
diff --git a/include/systems/implicit_system.h b/include/systems/implicit_system.h
index 8eef486c107..4b2e533316f 100644
--- a/include/systems/implicit_system.h
+++ b/include/systems/implicit_system.h
@@ -141,6 +141,7 @@ class ImplicitSystem : public ExplicitSystem
   virtual std::pair<unsigned int, Real>
   get_linear_solve_parameters() const;
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
   /**
    * Currently a no-op.
    *
@@ -149,6 +150,7 @@ class ImplicitSystem : public ExplicitSystem
    * pointer.
    */
   virtual void release_linear_solver(LinearSolver<Number> *) const;
+#endif // LIBMESH_ENABLE_DEPRECATED
 
   /**
    * Assembles a residual in \p rhs and/or a jacobian in \p matrix,
diff --git a/src/apps/embedding.C b/src/apps/embedding.C
index 22c8f8c1128..38883da9981 100644
--- a/src/apps/embedding.C
+++ b/src/apps/embedding.C
@@ -41,7 +41,6 @@ using namespace libMesh;
 void usage_error(const char * progname)
 {
   libMesh::out << "Options: " << progname << '\n'
-               << " --dim d               elem dimension\n"
                << " --elem type           elem type (e.g. TET14)\n"
                << " --childnum num        child number\n"
                << " --denominator num     denominator to use\n"
@@ -74,9 +73,6 @@ int main(int argc, char ** argv)
 
   GetPot cl(argc, argv);
 
-  const int dim =
-    assert_argument(cl, "--dim", argv[0], 0);
-
   const std::string elem_type_string =
     assert_argument(cl, "--elem", argv[0], std::string(""));
 
@@ -162,7 +158,7 @@ int main(int argc, char ** argv)
       std::cout << "      {";
       for (auto j : make_range(n_nodes))
         {
-          Real shape = FEInterface::shape(dim, fe_type, elem_type, j, pt);
+          Real shape = FEInterface::shape(fe_type, elem.get(), j, pt);
 
           // Don't print -0 or 1e-17; we don't tolerate FP error at 0
           if (std::abs(shape) < TOLERANCE*std::sqrt(TOLERANCE))
diff --git a/src/fe/fe.C b/src/fe/fe.C
index 82069be3f33..ef48b4eca28 100644
--- a/src/fe/fe.C
+++ b/src/fe/fe.C
@@ -60,7 +60,9 @@ FE<Dim,T>::FE (const FEType & fet) :
 template <unsigned int Dim, FEFamily T>
 unsigned int FE<Dim,T>::n_shape_functions () const
 {
-  return FE<Dim,T>::n_dofs (this->elem_type,
+  libmesh_deprecated(); // Use n_dofs(Elem *)
+
+  return FE<Dim,T>::n_dofs (this->get_type(),
                             this->fe_type.order + this->_p_level);
 }
 
@@ -100,7 +102,7 @@ void FE<Dim,T>::dofs_on_side(const Elem * const elem,
   for (unsigned int n = 0; n != n_nodes; ++n)
     {
       const unsigned int n_dofs =
-          n_dofs_at_node(elem->type(), static_cast<Order>(o + add_p_level*elem->p_level()), n);
+          n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
       if (elem->is_node_on_side(n, s))
         for (unsigned int i = 0; i != n_dofs; ++i)
           di.push_back(nodenum++);
@@ -127,7 +129,7 @@ void FE<Dim,T>::dofs_on_edge(const Elem * const elem,
   for (unsigned int n = 0; n != n_nodes; ++n)
     {
       const unsigned int n_dofs =
-          n_dofs_at_node(elem->type(), static_cast<Order>(o + add_p_level*elem->p_level()), n);
+          n_dofs_at_node(*elem, static_cast<Order>(o + add_p_level*elem->p_level()), n);
       if (elem->is_node_on_edge(n, e))
         for (unsigned int i = 0; i != n_dofs; ++i)
           di.push_back(nodenum++);
@@ -333,8 +335,7 @@ void FE<Dim,T>::reinit_dual_shape_coeffs(const Elem * elem,
   this->_p_level = this->_add_p_level_in_reinit * elem->p_level();
 
   const unsigned int n_shapes =
-    this->n_shape_functions(this->get_type(),
-                            this->get_order());
+    this->n_dofs(elem, this->get_order());
 
   std::vector<std::vector<OutputShape>> phi_vals;
   phi_vals.resize(n_shapes);
@@ -405,8 +406,7 @@ void FE<Dim,T>::init_shape_functions(const std::vector<Point> & qp,
   // Number of shape functions in the finite element approximation
   // space.
   const unsigned int n_approx_shape_functions =
-    this->n_shape_functions(this->get_type(),
-                            this->get_order());
+    this->n_dofs(elem, this->get_order());
 
   // Maybe we already have correctly-sized data?  Check data sizes,
   // and get ready to break out of a "loop" if all these resize()
diff --git a/src/fe/fe_abstract.C b/src/fe/fe_abstract.C
index a56b7ead6e6..bea2ac4d78a 100644
--- a/src/fe/fe_abstract.C
+++ b/src/fe/fe_abstract.C
@@ -623,8 +623,14 @@ void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point>
     }
 }
 
+
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
 bool FEAbstract::on_reference_element(const Point & p, const ElemType t, const Real eps)
 {
+  // Use Elem::on_reference_element() instead
+  libmesh_deprecated();
+
   libmesh_assert_greater_equal (eps, 0.);
 
   const Real xi   = p(0);
@@ -825,6 +831,7 @@ bool FEAbstract::on_reference_element(const Point & p, const ElemType t, const R
 
   return false;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/fe/fe_bernstein.C b/src/fe/fe_bernstein.C
index 86cc6de2beb..1d79ae09f2a 100644
--- a/src/fe/fe_bernstein.C
+++ b/src/fe/fe_bernstein.C
@@ -109,7 +109,7 @@ void bernstein_nodal_soln(const Elem * elem,
 
 
 
-unsigned int bernstein_n_dofs(const ElemType t, const Order o)
+unsigned int BERNSTEIN_n_dofs(const ElemType t, const Order o)
 {
   switch (t)
     {
@@ -168,12 +168,19 @@ unsigned int bernstein_n_dofs(const ElemType t, const Order o)
     default:
       libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
     }
-} // bernstein_n_dofs()
+} // BERNSTEIN_n_dofs()
 
 
 
+unsigned int BERNSTEIN_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return BERNSTEIN_n_dofs(e->type(), o);
+}
+
 
-unsigned int bernstein_n_dofs_at_node(const ElemType t,
+
+unsigned int BERNSTEIN_n_dofs_at_node(const ElemType t,
                                       const Order o,
                                       const unsigned int n)
 {
@@ -340,12 +347,20 @@ unsigned int bernstein_n_dofs_at_node(const ElemType t,
     default:
       libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
     }
-} // bernstein_n_dofs_at_node()
+} // BERNSTEIN_n_dofs_at_node()
+
 
 
+unsigned int BERNSTEIN_n_dofs_at_node(const Elem & e,
+                                      const Order o,
+                                      const unsigned int n)
+{
+  return BERNSTEIN_n_dofs_at_node(e.type(), o, n);
+}
+
 
 
-unsigned int bernstein_n_dofs_per_elem(const ElemType t, const Order o)
+unsigned int BERNSTEIN_n_dofs_per_elem(const ElemType t, const Order o)
 {
   switch (t)
     {
@@ -390,7 +405,14 @@ unsigned int bernstein_n_dofs_per_elem(const ElemType t, const Order o)
     default:
       libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
     }
-} // bernstein_n_dofs_per_elem
+} // BERNSTEIN_n_dofs_per_elem
+
+
+
+unsigned int BERNSTEIN_n_dofs_per_elem(const Elem & e, const Order o)
+{
+  return BERNSTEIN_n_dofs_per_elem(e.type(), o);
+}
 
 } // anonymous namespace
 
@@ -399,24 +421,8 @@ unsigned int bernstein_n_dofs_per_elem(const ElemType t, const Order o)
 LIBMESH_FE_NODAL_SOLN(BERNSTEIN, bernstein_nodal_soln)
 LIBMESH_FE_SIDE_NODAL_SOLN(BERNSTEIN)
 
-
-// Full specialization of n_dofs() function for every dimension
-template <> unsigned int FE<0,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
-template <> unsigned int FE<1,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
-template <> unsigned int FE<2,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
-template <> unsigned int FE<3,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
-
-// Full specialization of n_dofs_at_node() function for every dimension.
-template <> unsigned int FE<0,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<1,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<2,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<3,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
-
-// Full specialization of n_dofs_per_elem() function for every dimension.
-template <> unsigned int FE<0,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<1,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<2,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<3,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
+// Instantiate n_dofs*() functions for every dimension
+LIBMESH_DEFAULT_NDOFS(BERNSTEIN)
 
 // Bernstein FEMs are C^0 continuous
 template <> FEContinuity FE<0,BERNSTEIN>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_clough.C b/src/fe/fe_clough.C
index f02412fc439..f273cde67d4 100644
--- a/src/fe/fe_clough.C
+++ b/src/fe/fe_clough.C
@@ -87,8 +87,7 @@ void clough_nodal_soln(const Elem * elem,
 
 
 
-
-unsigned int clough_n_dofs(const ElemType t, const Order o)
+unsigned int CLOUGH_n_dofs(const ElemType t, const Order o)
 {
   if (t == INVALID_ELEM)
     return 0;
@@ -125,13 +124,19 @@ unsigned int clough_n_dofs(const ElemType t, const Order o)
     default:
       libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
     }
-} // clough_n_dofs()
+} // CLOUGH_n_dofs()
+
 
 
+unsigned int CLOUGH_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return CLOUGH_n_dofs(e->type(), o);
+}
 
 
 
-unsigned int clough_n_dofs_at_node(const ElemType t,
+unsigned int CLOUGH_n_dofs_at_node(const ElemType t,
                                    const Order o,
                                    const unsigned int n)
 {
@@ -204,12 +209,20 @@ unsigned int clough_n_dofs_at_node(const ElemType t,
     default:
       libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
     }
-} // clough_n_dofs_at_node()
+} // CLOUGH_n_dofs_at_node()
+
 
 
+unsigned int CLOUGH_n_dofs_at_node(const Elem & e,
+                                   const Order o,
+                                   const unsigned int n)
+{
+  return CLOUGH_n_dofs_at_node(e.type(), o, n);
+}
+
 
 
-unsigned int clough_n_dofs_per_elem(const ElemType t, const Order o)
+unsigned int CLOUGH_n_dofs_per_elem(const ElemType t, const Order o)
 {
   switch (o)
     {
@@ -233,7 +246,14 @@ unsigned int clough_n_dofs_per_elem(const ElemType t, const Order o)
     default:
       libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
     }
-} // clough_n_dofs_per_elem()
+} // CLOUGH_n_dofs_per_elem()
+
+
+
+unsigned int CLOUGH_n_dofs_per_elem(const Elem & e, const Order o)
+{
+  return CLOUGH_n_dofs_per_elem(e.type(), o);
+}
 
 } // anonymous
 
@@ -245,24 +265,9 @@ LIBMESH_FE_NODAL_SOLN(CLOUGH, clough_nodal_soln)
 LIBMESH_FE_SIDE_NODAL_SOLN(CLOUGH)
 
 
-// Full specialization of n_dofs() function for every dimension
-template <> unsigned int FE<0,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
-template <> unsigned int FE<1,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
-template <> unsigned int FE<2,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
-template <> unsigned int FE<3,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
-
-
-// Full specialization of n_dofs_at_node() function for every dimension.
-template <> unsigned int FE<0,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<1,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<2,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<3,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
+// Instantiate n_dofs*() functions for every dimension
+LIBMESH_DEFAULT_NDOFS(CLOUGH)
 
-// Full specialization of n_dofs_per_elem() function for every dimension.
-template <> unsigned int FE<0,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<1,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<2,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<3,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
 
 // Clough FEMs are C^1 continuous
 template <> FEContinuity FE<0,CLOUGH>::get_continuity() const { return C_ONE; }
diff --git a/src/fe/fe_hermite.C b/src/fe/fe_hermite.C
index 7c9a5c458a8..fc36714852a 100644
--- a/src/fe/fe_hermite.C
+++ b/src/fe/fe_hermite.C
@@ -71,7 +71,7 @@ void hermite_nodal_soln(const Elem * elem,
 
 
 
-unsigned int hermite_n_dofs(const ElemType t, const Order o)
+unsigned int HERMITE_n_dofs(const ElemType t, const Order o)
 {
 #ifdef DEBUG
   libmesh_error_msg_if(o < 3, "Error: Hermite elements require order>=3, but you asked for order=" << o);
@@ -111,12 +111,19 @@ unsigned int hermite_n_dofs(const ElemType t, const Order o)
     default:
       libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
     }
-} // hermite_n_dofs()
+} // HERMITE_n_dofs()
 
 
 
+unsigned int HERMITE_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return HERMITE_n_dofs(e->type(), o);
+}
+
 
-unsigned int hermite_n_dofs_at_node(const ElemType t,
+
+unsigned int HERMITE_n_dofs_at_node(const ElemType t,
                                     const Order o,
                                     const unsigned int n)
 {
@@ -233,11 +240,20 @@ unsigned int hermite_n_dofs_at_node(const ElemType t,
     default:
       libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
     }
-} // hermite_n_dofs_at_node()
+} // HERMITE_n_dofs_at_node()
+
 
 
+unsigned int HERMITE_n_dofs_at_node(const Elem & e,
+                                    const Order o,
+                                    const unsigned int n)
+{
+  return HERMITE_n_dofs_at_node(e.type(), o, n);
+}
+
 
-unsigned int hermite_n_dofs_per_elem(const ElemType t,
+
+unsigned int HERMITE_n_dofs_per_elem(const ElemType t,
                                      const Order o)
 {
   libmesh_assert_greater (o, 2);
@@ -271,7 +287,15 @@ unsigned int hermite_n_dofs_per_elem(const ElemType t,
     default:
       libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
     }
-} // hermite_n_dofs_per_elem()
+} // HERMITE_n_dofs_per_elem()
+
+
+
+unsigned int HERMITE_n_dofs_per_elem(const Elem & e,
+                                     const Order o)
+{
+  return HERMITE_n_dofs_per_elem(e.type(), o);
+}
 
 
 } // anonymous namespace
@@ -282,26 +306,9 @@ LIBMESH_FE_NODAL_SOLN(HERMITE, hermite_nodal_soln)
 LIBMESH_FE_SIDE_NODAL_SOLN(HERMITE)
 
 
-// Do full-specialization for every dimension, instead
-// of explicit instantiation at the end of this function.
-// This could be macro-ified.
-template <> unsigned int FE<0,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
-template <> unsigned int FE<1,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
-template <> unsigned int FE<2,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
-template <> unsigned int FE<3,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
-
-// Do full-specialization for every dimension, instead
-// of explicit instantiation at the end of this function.
-template <> unsigned int FE<0,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<1,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<2,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<3,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
-
-// Full specialization of n_dofs_per_elem() function for every dimension.
-template <> unsigned int FE<0,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<1,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<2,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<3,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
+// Instantiate n_dofs*() functions for every dimension
+LIBMESH_DEFAULT_NDOFS(HERMITE)
+
 
 // Hermite FEMs are C^1 continuous
 template <> FEContinuity FE<0,HERMITE>::get_continuity() const { return C_ONE; }
diff --git a/src/fe/fe_hierarchic.C b/src/fe/fe_hierarchic.C
index 856c1b4d5c8..288f9dc06c7 100644
--- a/src/fe/fe_hierarchic.C
+++ b/src/fe/fe_hierarchic.C
@@ -97,8 +97,7 @@ void hierarchic_nodal_soln(const Elem * elem,
 
 
 
-
-unsigned int hierarchic_n_dofs(const ElemType t, const Order o)
+unsigned int HIERARCHIC_n_dofs(const ElemType t, const Order o)
 {
   libmesh_assert_greater (o, 0);
   switch (t)
@@ -152,12 +151,19 @@ unsigned int hierarchic_n_dofs(const ElemType t, const Order o)
     default:
       libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for HIERARCHIC FE family!");
     }
-} // hierarchic_n_dofs()
+} // HIERARCHIC_n_dofs()
+
 
 
+unsigned int HIERARCHIC_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return HIERARCHIC_n_dofs(e->type(), o);
+}
+
 
 
-unsigned int hierarchic_n_dofs_at_node(const ElemType t,
+unsigned int HIERARCHIC_n_dofs_at_node(const ElemType t,
                                        const Order o,
                                        const unsigned int n)
 {
@@ -382,12 +388,20 @@ unsigned int hierarchic_n_dofs_at_node(const ElemType t,
     default:
       libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
     }
-} // hierarchic_n_dofs_at_node()
+} // HIERARCHIC_n_dofs_at_node()
 
 
 
+unsigned int HIERARCHIC_n_dofs_at_node(const Elem & e,
+                                       const Order o,
+                                       const unsigned int n)
+{
+  return HIERARCHIC_n_dofs_at_node(e.type(), o, n);
+}
+
 
-unsigned int hierarchic_n_dofs_per_elem(const ElemType t,
+
+unsigned int HIERARCHIC_n_dofs_per_elem(const ElemType t,
                                         const Order o)
 {
   libmesh_assert_greater (o, 0);
@@ -442,7 +456,15 @@ unsigned int hierarchic_n_dofs_per_elem(const ElemType t,
     default:
       libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
     }
-} // hierarchic_n_dofs_per_elem()
+} // HIERARCHIC_n_dofs_per_elem()
+
+
+
+unsigned int HIERARCHIC_n_dofs_per_elem(const Elem & e,
+                                        const Order o)
+{
+  return HIERARCHIC_n_dofs_per_elem(e.type(), o);
+}
 
 } // anonymous namespace
 
@@ -452,23 +474,9 @@ LIBMESH_FE_NODAL_SOLN(HIERARCHIC, hierarchic_nodal_soln)
 LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC)
 
 
-// Full specialization of n_dofs() function for every dimension
-template <> unsigned int FE<0,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
-template <> unsigned int FE<1,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
-template <> unsigned int FE<2,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
-template <> unsigned int FE<3,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
-
-// Full specialization of n_dofs_at_node() function for every dimension.
-template <> unsigned int FE<0,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<1,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<2,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
-template <> unsigned int FE<3,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
+// Instantiate n_dofs*() functions for every dimension
+LIBMESH_DEFAULT_NDOFS(HIERARCHIC)
 
-// Full specialization of n_dofs_per_elem() function for every dimension.
-template <> unsigned int FE<0,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<1,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<2,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
-template <> unsigned int FE<3,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
 
 // Hierarchic FEMs are C^0 continuous
 template <> FEContinuity FE<0,HIERARCHIC>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_hierarchic_vec.C b/src/fe/fe_hierarchic_vec.C
index e44de992e4b..eceee989909 100644
--- a/src/fe/fe_hierarchic_vec.C
+++ b/src/fe/fe_hierarchic_vec.C
@@ -773,43 +773,42 @@ template <> RealGradient FE<3,L2_HIERARCHIC_VEC>::shape_second_deriv(const Elem
 
 #endif
 
-// Do full-specialization for every dimension, instead
-// of explicit instantiation at the end of this function.
-// This could be macro-ified.
-template <> unsigned int FE<0,HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,HIERARCHIC>::n_dofs(t,o); }
-template <> unsigned int FE<1,HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,HIERARCHIC>::n_dofs(t,o); }
-template <> unsigned int FE<2,HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,HIERARCHIC>::n_dofs(t,o); }
-template <> unsigned int FE<3,HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,HIERARCHIC>::n_dofs(t,o); }
 
+// Instantiate n_dofs*() functions for every dimension
+LIBMESH_DEFAULT_VEC_NDOFS(HIERARCHIC)
+
+// L2 Hierarchic has a simpler but non-default implementation, with
+// all their dofs on the element
 template <> unsigned int FE<0,L2_HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,L2_HIERARCHIC>::n_dofs(t,o); }
 template <> unsigned int FE<1,L2_HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,L2_HIERARCHIC>::n_dofs(t,o); }
 template <> unsigned int FE<2,L2_HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,L2_HIERARCHIC>::n_dofs(t,o); }
 template <> unsigned int FE<3,L2_HIERARCHIC_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,L2_HIERARCHIC>::n_dofs(t,o); }
 
-
-// Do full-specialization for every dimension, instead
-// of explicit instantiation at the end of this function.
-template <> unsigned int FE<0,HIERARCHIC_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,HIERARCHIC>::n_dofs_at_node(t,o,n); }
-template <> unsigned int FE<1,HIERARCHIC_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,HIERARCHIC>::n_dofs_at_node(t,o,n); }
-template <> unsigned int FE<2,HIERARCHIC_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 2*FE<2,HIERARCHIC>::n_dofs_at_node(t,o,n); }
-template <> unsigned int FE<3,HIERARCHIC_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<3,HIERARCHIC>::n_dofs_at_node(t,o,n); }
+template <> unsigned int FE<0,L2_HIERARCHIC_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,L2_HIERARCHIC>::n_dofs(e,o); }
+template <> unsigned int FE<1,L2_HIERARCHIC_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,L2_HIERARCHIC>::n_dofs(e,o); }
+template <> unsigned int FE<2,L2_HIERARCHIC_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,L2_HIERARCHIC>::n_dofs(e,o); }
+template <> unsigned int FE<3,L2_HIERARCHIC_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,L2_HIERARCHIC>::n_dofs(e,o); }
 
 template <> unsigned int FE<0,L2_HIERARCHIC_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<1,L2_HIERARCHIC_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<2,L2_HIERARCHIC_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,L2_HIERARCHIC_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
-template <> unsigned int FE<0,HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,HIERARCHIC>::n_dofs_per_elem(t,o); }
-template <> unsigned int FE<1,HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,HIERARCHIC>::n_dofs_per_elem(t,o); }
-template <> unsigned int FE<2,HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 2*FE<2,HIERARCHIC>::n_dofs_per_elem(t,o); }
-template <> unsigned int FE<3,HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 3*FE<3,HIERARCHIC>::n_dofs_per_elem(t,o); }
+template <> unsigned int FE<0,L2_HIERARCHIC_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,L2_HIERARCHIC_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,L2_HIERARCHIC_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,L2_HIERARCHIC_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
 
-// L2 Hierarchic elements have all their dofs on the element
 template <> unsigned int FE<0,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,L2_HIERARCHIC_VEC>::n_dofs(t, o); }
 template <> unsigned int FE<1,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,L2_HIERARCHIC_VEC>::n_dofs(t, o); }
 template <> unsigned int FE<2,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<2,L2_HIERARCHIC_VEC>::n_dofs(t, o); }
 template <> unsigned int FE<3,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<3,L2_HIERARCHIC_VEC>::n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,L2_HIERARCHIC_VEC>::n_dofs(&e, o); }
+template <> unsigned int FE<1,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,L2_HIERARCHIC_VEC>::n_dofs(&e, o); }
+template <> unsigned int FE<2,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<2,L2_HIERARCHIC_VEC>::n_dofs(&e, o); }
+template <> unsigned int FE<3,L2_HIERARCHIC_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<3,L2_HIERARCHIC_VEC>::n_dofs(&e, o); }
+
 // Hierarchic FEMs are always C^0 continuous
 template <> FEContinuity FE<0,HIERARCHIC_VEC>::get_continuity() const { return C_ZERO; }
 template <> FEContinuity FE<1,HIERARCHIC_VEC>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_interface.C b/src/fe/fe_interface.C
index 8bfe98e8298..4e5caa9b98c 100644
--- a/src/fe/fe_interface.C
+++ b/src/fe/fe_interface.C
@@ -268,12 +268,12 @@ FEInterface::is_InfFE_elem(const ElemType et)
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::n_shape_functions(const unsigned int dim,
                                             const FEType & fe_t,
                                             const ElemType t)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
   /*
@@ -291,6 +291,7 @@ unsigned int FEInterface::n_shape_functions(const unsigned int dim,
 
   fe_with_vec_switch(n_shape_functions(t, o));
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -317,7 +318,7 @@ FEInterface::n_shape_functions(const FEType & fe_t,
   // Account for Elem::p_level() when computing total_order
   auto total_order = fe_t.order + add_p_level*elem->p_level();
 
-  fe_with_vec_switch(n_shape_functions(elem->type(), total_order));
+  fe_with_vec_switch(n_dofs(elem, total_order));
 }
 
 
@@ -345,11 +346,12 @@ FEInterface::n_shape_functions(const FEType & fe_t,
   // Ignore Elem::p_level() and instead use extra_order to compute total_order.
   auto total_order = fe_t.order + extra_order;
 
-  fe_with_vec_switch(n_shape_functions(elem->type(), total_order));
+  fe_with_vec_switch(n_dofs(elem, total_order));
 }
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::n_dofs(const unsigned int dim,
                                  const FEType & fe_t,
                                  const ElemType t)
@@ -378,10 +380,9 @@ FEInterface::n_dofs (const unsigned int dim,
 {
   libmesh_deprecated();
 
-  FEType p_refined_fe_t = fe_t;
-  p_refined_fe_t.order = p_refined_fe_t.order + elem->p_level();
-  return FEInterface::n_dofs(dim, p_refined_fe_t, elem->type());
+  fe_with_vec_switch(n_dofs(elem, fe_t.order + elem->p_level()));
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -402,9 +403,7 @@ FEInterface::n_dofs(const FEType & fe_t,
 #endif
 
   // Account for Elem::p_level() when computing total_order
-  auto total_order = fe_t.order + add_p_level*elem->p_level();
-
-  fe_with_vec_switch(n_dofs(elem->type(), total_order));
+  fe_with_vec_switch(n_dofs(elem, fe_t.order + add_p_level*elem->p_level()));
 }
 
 
@@ -428,18 +427,18 @@ FEInterface::n_dofs(const FEType & fe_t,
   // Elem::p_level() is ignored, extra_order is used instead.
   auto total_order = fe_t.order + extra_order;
 
-  fe_with_vec_switch(n_dofs(elem->type(), total_order));
+  fe_with_vec_switch(n_dofs(elem, total_order));
 }
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::n_dofs_at_node(const unsigned int dim,
                                          const FEType & fe_t,
                                          const ElemType t,
                                          const unsigned int n)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -463,6 +462,7 @@ FEInterface::n_dofs_at_node_function(const unsigned int dim,
 
   fe_with_vec_switch(n_dofs_at_node);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -497,7 +497,7 @@ FEInterface::n_dofs_at_node(const FEType & fe_t,
   // Account for Elem::p_level() when computing total_order
   auto total_order = fe_t.order + add_p_level*elem->p_level();
 
-  fe_with_vec_switch(n_dofs_at_node(elem->type(), total_order, n));
+  fe_with_vec_switch(n_dofs_at_node(*elem, total_order, n));
 }
 
 
@@ -521,17 +521,17 @@ FEInterface::n_dofs_at_node(const FEType & fe_t,
   // Ignore Elem::p_level() and instead use extra_order to compute total_order.
   auto total_order = fe_t.order + extra_order;
 
-  fe_with_vec_switch(n_dofs_at_node(elem->type(), total_order, n));
+  fe_with_vec_switch(n_dofs_at_node(*elem, total_order, n));
 }
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim,
                                           const FEType & fe_t,
                                           const ElemType t)
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -544,6 +544,7 @@ unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim,
 
   fe_with_vec_switch(n_dofs_per_elem(t, o));
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -565,7 +566,7 @@ FEInterface::n_dofs_per_elem(const FEType & fe_t,
   // Account for Elem::p_level() when computing total_order
   auto total_order = fe_t.order + add_p_level*elem->p_level();
 
-  fe_with_vec_switch(n_dofs_per_elem(elem->type(), total_order));
+  fe_with_vec_switch(n_dofs_per_elem(*elem, total_order));
 }
 
 
@@ -588,7 +589,7 @@ FEInterface::n_dofs_per_elem(const FEType & fe_t,
   // Ignore Elem::p_level() and instead use extra_order to compute total_order.
   auto total_order = fe_t.order + extra_order;
 
-  fe_with_vec_switch(n_dofs_per_elem(elem->type(), total_order));
+  fe_with_vec_switch(n_dofs_per_elem(*elem, total_order));
 }
 
 
@@ -671,11 +672,13 @@ void FEInterface::side_nodal_soln(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 Point FEInterface::map(unsigned int dim,
                        const FEType & fe_t,
                        const Elem * elem,
                        const Point & p)
 {
+  libmesh_deprecated();
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
   if (is_InfFE_elem(elem->type()))
     return ifem_map(dim, fe_t, elem, p);
@@ -685,8 +688,6 @@ Point FEInterface::map(unsigned int dim,
 
 
 
-
-
 Point FEInterface::inverse_map (const unsigned int dim,
                                 const FEType & fe_t,
                                 const Elem * elem,
@@ -694,6 +695,7 @@ Point FEInterface::inverse_map (const unsigned int dim,
                                 const Real tolerance,
                                 const bool secure)
 {
+  libmesh_deprecated();
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
   if (is_InfFE_elem(elem->type()))
@@ -706,7 +708,6 @@ Point FEInterface::inverse_map (const unsigned int dim,
 
 
 
-
 void FEInterface::inverse_map (const unsigned int dim,
                                const FEType & fe_t,
                                const Elem * elem,
@@ -715,6 +716,8 @@ void FEInterface::inverse_map (const unsigned int dim,
                                const Real tolerance,
                                const bool secure)
 {
+  libmesh_deprecated();
+
   const std::size_t n_pts = physical_points.size();
 
   // Resize the vector
@@ -753,15 +756,13 @@ bool FEInterface::on_reference_element(const Point & p,
 
 
 
-
 Real FEInterface::shape(const unsigned int dim,
                         const FEType & fe_t,
                         const ElemType t,
                         const unsigned int i,
                         const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -781,8 +782,7 @@ Real FEInterface::shape(const unsigned int dim,
                         const unsigned int i,
                         const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -795,6 +795,7 @@ Real FEInterface::shape(const unsigned int dim,
 
   fe_switch(shape(elem,o,i,p));
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -856,6 +857,7 @@ FEInterface::shape(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template<>
 void FEInterface::shape<Real>(const unsigned int dim,
                               const FEType & fe_t,
@@ -864,8 +866,7 @@ void FEInterface::shape<Real>(const unsigned int dim,
                               const Point & p,
                               Real & phi)
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -900,6 +901,8 @@ void FEInterface::shape<Real>(const unsigned int dim,
   return;
 }
 
+
+
 template<>
 void FEInterface::shape<Real>(const unsigned int dim,
                               const FEType & fe_t,
@@ -908,8 +911,7 @@ void FEInterface::shape<Real>(const unsigned int dim,
                               const Point & p,
                               Real & phi)
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 
@@ -942,6 +944,7 @@ void FEInterface::shape<Real>(const unsigned int dim,
 
   return;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1127,8 +1130,7 @@ void FEInterface::all_shapes<Real>(const unsigned int dim,
 
 
 
-
-
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template<>
 void FEInterface::shape<RealGradient>(const unsigned int dim,
                                       const FEType & fe_t,
@@ -1137,8 +1139,7 @@ void FEInterface::shape<RealGradient>(const unsigned int dim,
                                       const Point & p,
                                       RealGradient & phi)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   // This API does not currently support infinite elements.
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
@@ -1167,6 +1168,7 @@ void FEInterface::shape<RealGradient>(const unsigned int dim,
       libmesh_error_msg("Invalid dimension = " << dim);
     }
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1337,13 +1339,13 @@ void FEInterface::all_shapes<RealGradient>(const unsigned int dim,
 }
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 FEInterface::shape_ptr
 FEInterface::shape_function(const unsigned int dim,
                             const FEType & fe_t,
                             const ElemType libmesh_inf_var(t))
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
   if (is_InfFE_elem(t))
@@ -1353,6 +1355,7 @@ FEInterface::shape_function(const unsigned int dim,
 #endif
   fe_switch(shape);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1371,6 +1374,8 @@ FEInterface::shape_function(const FEType & fe_t,
 }
 
 
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
 Real FEInterface::shape_deriv(const unsigned int dim,
                               const FEType & fe_t,
                               const ElemType t,
@@ -1378,8 +1383,7 @@ Real FEInterface::shape_deriv(const unsigned int dim,
                               const unsigned int j,
                               const Point & p)
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   libmesh_assert_greater (dim,j);
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
@@ -1403,8 +1407,7 @@ Real FEInterface::shape_deriv(const unsigned int dim,
                               const unsigned int j,
                               const Point & p)
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   libmesh_assert_greater (dim,j);
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
@@ -1432,6 +1435,7 @@ Real FEInterface::shape_deriv(const unsigned int dim,
     }
   return 0;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1621,13 +1625,13 @@ void FEInterface::shape_derivs<RealGradient>(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 FEInterface::shape_deriv_ptr
 FEInterface::shape_deriv_function(const unsigned int dim,
                                   const FEType & fe_t,
                                   const ElemType libmesh_inf_var(t))
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
   if (is_InfFE_elem(t))
@@ -1637,6 +1641,7 @@ FEInterface::shape_deriv_function(const unsigned int dim,
 #endif
   fe_switch(shape_deriv);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1658,6 +1663,7 @@ FEInterface::shape_deriv_function(const FEType & fe_t,
 
 
 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
+#ifdef LIBMESH_ENABLE_DEPRECATED
 Real FEInterface::shape_second_deriv(const unsigned int dim,
                                      const FEType & fe_t,
                                      const ElemType t,
@@ -1665,8 +1671,7 @@ Real FEInterface::shape_second_deriv(const unsigned int dim,
                                      const unsigned int j,
                                      const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   libmesh_assert_greater_equal (dim*(dim-1),j);
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
@@ -1700,8 +1705,7 @@ Real FEInterface::shape_second_deriv(const unsigned int dim,
                                      const unsigned int j,
                                      const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   libmesh_assert_greater_equal (dim*(dim-1),j);
 
@@ -1725,6 +1729,7 @@ Real FEInterface::shape_second_deriv(const unsigned int dim,
     }
   return 0;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1806,13 +1811,13 @@ Real FEInterface::shape_second_deriv(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 FEInterface::shape_second_deriv_ptr
 FEInterface::shape_second_deriv_function(const unsigned int dim,
                                          const FEType & fe_t,
                                          const ElemType libmesh_inf_var(t))
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
   if (is_InfFE_elem(t))
@@ -1820,6 +1825,7 @@ FEInterface::shape_second_deriv_function(const unsigned int dim,
 #endif
   fe_switch(shape_second_deriv);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -1839,6 +1845,7 @@ FEInterface::shape_second_deriv_function(const FEType & fe_t,
 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template<>
 void FEInterface::shape<RealGradient>(const unsigned int dim,
                                       const FEType & fe_t,
@@ -1847,8 +1854,7 @@ void FEInterface::shape<RealGradient>(const unsigned int dim,
                                       const Point & p,
                                       RealGradient & phi)
 {
-  // TODO
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   // This is actually an issue for infinite elements: They require type 'Gradient'!
   if (elem->infinite())
@@ -1876,6 +1882,8 @@ void FEInterface::shape<RealGradient>(const unsigned int dim,
 
   return;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
+
 
 void FEInterface::compute_data(const unsigned int dim,
                                const FEType & fe_t,
@@ -1938,11 +1946,11 @@ void FEInterface::compute_data(const unsigned int dim,
       // shape_deriv() functions since they have access to the elem
       // pointer. Note that we are already using the n_dof value
       // appropriate to the elevated p-level.
-      data.shape[n] = shape(dim, fe_t, elem, n, p);
+      data.shape[n] = shape(fe_t, elem, n, p);
       if (data.need_derivative())
         {
           for (unsigned int j=0; j<dim; j++)
-            data.dshape[n](j) = shape_deriv(dim, fe_t, elem, n, j, p);
+            data.dshape[n](j) = shape_deriv(fe_t, elem, n, j, p);
         }
     }
 }
diff --git a/src/fe/fe_interface_inf_fe.C b/src/fe/fe_interface_inf_fe.C
index ecec5dca115..3046ee0a9e1 100644
--- a/src/fe/fe_interface_inf_fe.C
+++ b/src/fe/fe_interface_inf_fe.C
@@ -36,12 +36,12 @@ namespace libMesh
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::ifem_n_shape_functions(const unsigned int dim,
                                                  const FEType & fe_t,
                                                  const ElemType t)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   switch (dim)
     {
@@ -66,6 +66,7 @@ unsigned int FEInterface::ifem_n_shape_functions(const unsigned int dim,
       libmesh_error_msg("Unsupported dim = " << dim);
     }
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -98,6 +99,7 @@ unsigned int FEInterface::ifem_n_shape_functions(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::ifem_n_dofs(const unsigned int dim,
                                       const FEType & fe_t,
                                       const ElemType t)
@@ -127,6 +129,7 @@ unsigned int FEInterface::ifem_n_dofs(const unsigned int dim,
       libmesh_error_msg("Unsupported dim = " << dim);
     }
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -160,13 +163,13 @@ FEInterface::ifem_n_dofs(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::ifem_n_dofs_at_node(const unsigned int dim,
                                               const FEType & fe_t,
                                               const ElemType t,
                                               const unsigned int n)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   switch (dim)
     {
@@ -191,6 +194,7 @@ unsigned int FEInterface::ifem_n_dofs_at_node(const unsigned int dim,
       libmesh_error_msg("Unsupported dim = " << dim);
     }
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -226,12 +230,12 @@ unsigned int FEInterface::ifem_n_dofs_at_node(const FEType & fe_t,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int FEInterface::ifem_n_dofs_per_elem(const unsigned int dim,
                                                const FEType & fe_t,
                                                const ElemType t)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   switch (dim)
     {
@@ -256,6 +260,7 @@ unsigned int FEInterface::ifem_n_dofs_per_elem(const unsigned int dim,
       libmesh_error_msg("Unsupported dim = " << dim);
     }
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -724,7 +729,7 @@ void FEInterface::ifem_inverse_map (const unsigned int dim,
 
 
 
-
+#ifdef LIBMESH_ENABLE_DEPRECATED
 bool FEInterface::ifem_on_reference_element(const Point & p,
                                             const ElemType t,
                                             const Real eps)
@@ -740,26 +745,26 @@ Real FEInterface::ifem_shape(const unsigned int dim,
                              const unsigned int i,
                              const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   inf_fe_switch(shape(fe_t, t, i, p));
 }
 
 
 
-
 Real FEInterface::ifem_shape(const unsigned int dim,
                              const FEType & fe_t,
                              const Elem * elem,
                              const unsigned int i,
                              const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   inf_fe_switch( shape(fe_t, elem, i, p));
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
+
+
 
 Real FEInterface::ifem_shape(const FEType & fe_t,
                              const Elem * elem,
@@ -772,6 +777,9 @@ Real FEInterface::ifem_shape(const FEType & fe_t,
   inf_fe_switch( shape(fe_t, elem, i, p));
 }
 
+
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
 Real FEInterface::ifem_shape_deriv (const unsigned int dim,
                                     const FEType & fe_t,
                                     const Elem * elem,
@@ -779,13 +787,13 @@ Real FEInterface::ifem_shape_deriv (const unsigned int dim,
                                     const unsigned int j,
                                     const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   inf_fe_switch(shape_deriv(fe_t, elem, i, j, p));
 }
 
 
+
 Real FEInterface::ifem_shape_deriv(const unsigned int dim,
                                    const FEType & fe_t,
                                    const ElemType t,
@@ -793,11 +801,11 @@ Real FEInterface::ifem_shape_deriv(const unsigned int dim,
                                    const unsigned int j,
                                    const Point & p)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   inf_fe_switch(shape_deriv(fe_t, t, i, j, p));
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/fe/fe_l2_hierarchic.C b/src/fe/fe_l2_hierarchic.C
index 8a9ead2c2d4..15bc758f495 100644
--- a/src/fe/fe_l2_hierarchic.C
+++ b/src/fe/fe_l2_hierarchic.C
@@ -97,7 +97,6 @@ void l2_hierarchic_nodal_soln(const Elem * elem,
 
 
 
-
 unsigned int l2_hierarchic_n_dofs(const ElemType t, const Order o)
 {
   libmesh_assert_greater (o, 0);
@@ -141,6 +140,14 @@ unsigned int l2_hierarchic_n_dofs(const ElemType t, const Order o)
 } // l2_hierarchic_n_dofs()
 
 
+
+unsigned int l2_hierarchic_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return l2_hierarchic_n_dofs(e->type(), o);
+}
+
+
 } // anonymous namespace
 
 
@@ -155,6 +162,11 @@ template <> unsigned int FE<1,L2_HIERARCHIC>::n_dofs(const ElemType t, const Ord
 template <> unsigned int FE<2,L2_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return l2_hierarchic_n_dofs(t, o); }
 template <> unsigned int FE<3,L2_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return l2_hierarchic_n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return l2_hierarchic_n_dofs(e, o); }
+template <> unsigned int FE<1,L2_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return l2_hierarchic_n_dofs(e, o); }
+template <> unsigned int FE<2,L2_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return l2_hierarchic_n_dofs(e, o); }
+template <> unsigned int FE<3,L2_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return l2_hierarchic_n_dofs(e, o); }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 // Discontinuous L2 elements only have interior nodes
 template <> unsigned int FE<0,L2_HIERARCHIC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
@@ -162,12 +174,22 @@ template <> unsigned int FE<1,L2_HIERARCHIC>::n_dofs_at_node(const ElemType, con
 template <> unsigned int FE<2,L2_HIERARCHIC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,L2_HIERARCHIC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,L2_HIERARCHIC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,L2_HIERARCHIC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,L2_HIERARCHIC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,L2_HIERARCHIC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 template <> unsigned int FE<0,L2_HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_hierarchic_n_dofs(t, o); }
 template <> unsigned int FE<1,L2_HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_hierarchic_n_dofs(t, o); }
 template <> unsigned int FE<2,L2_HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_hierarchic_n_dofs(t, o); }
 template <> unsigned int FE<3,L2_HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_hierarchic_n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_HIERARCHIC>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_hierarchic_n_dofs(&e, o); }
+template <> unsigned int FE<1,L2_HIERARCHIC>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_hierarchic_n_dofs(&e, o); }
+template <> unsigned int FE<2,L2_HIERARCHIC>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_hierarchic_n_dofs(&e, o); }
+template <> unsigned int FE<3,L2_HIERARCHIC>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_hierarchic_n_dofs(&e, o); }
+
 // L2 Hierarchic FEMs are C^0 continuous
 template <> FEContinuity FE<0,L2_HIERARCHIC>::get_continuity() const { return DISCONTINUOUS; }
 template <> FEContinuity FE<1,L2_HIERARCHIC>::get_continuity() const { return DISCONTINUOUS; }
diff --git a/src/fe/fe_l2_lagrange.C b/src/fe/fe_l2_lagrange.C
index 406f8ee63a8..a745d87e9fe 100644
--- a/src/fe/fe_l2_lagrange.C
+++ b/src/fe/fe_l2_lagrange.C
@@ -201,6 +201,14 @@ unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o)
     }
 }
 
+
+
+unsigned int l2_lagrange_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return l2_lagrange_n_dofs(e->type(), o);
+}
+
 } // anonymous namespace
 
 
@@ -217,6 +225,10 @@ template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs(const ElemType t, const Order
 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_n_dofs(e, o); }
+template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_n_dofs(e, o); }
+template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_n_dofs(e, o); }
+template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs(const Elem * e, const Order o) { return l2_lagrange_n_dofs(e, o); }
 
 // L2 Lagrange elements have all dofs on elements (hence none on nodes)
 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
@@ -224,6 +236,10 @@ template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const
 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
 
 // L2 Lagrange elements have all dofs on elements
 template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
@@ -231,6 +247,11 @@ template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, co
 template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
 template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_n_dofs(&e, o); }
+template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_n_dofs(&e, o); }
+template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_n_dofs(&e, o); }
+template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_per_elem(const Elem & e, const Order o) { return l2_lagrange_n_dofs(&e, o); }
+
 // L2 Lagrange FEMs are DISCONTINUOUS
 template <> FEContinuity FE<0,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; }
 template <> FEContinuity FE<1,L2_LAGRANGE>::get_continuity() const { return DISCONTINUOUS; }
diff --git a/src/fe/fe_lagrange.C b/src/fe/fe_lagrange.C
index 94e64b585a2..f35e621193b 100644
--- a/src/fe/fe_lagrange.C
+++ b/src/fe/fe_lagrange.C
@@ -603,6 +603,13 @@ unsigned int lagrange_n_dofs(const ElemType t, const Order o)
 
 
 
+unsigned int lagrange_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return lagrange_n_dofs(e->type(), o);
+}
+
+
 
 unsigned int lagrange_n_dofs_at_node(const ElemType t,
                                      const Order o,
@@ -841,6 +848,15 @@ unsigned int lagrange_n_dofs_at_node(const ElemType t,
 }
 
 
+
+unsigned int lagrange_n_dofs_at_node(const Elem & e,
+                                     const Order o,
+                                     const unsigned int n)
+{
+  return lagrange_n_dofs_at_node(e.type(), o, n);
+}
+
+
 #ifdef LIBMESH_ENABLE_AMR
 void lagrange_compute_constraints (DofConstraints & constraints,
                                    DofMap & dof_map,
@@ -1003,8 +1019,7 @@ void lagrange_compute_constraints (DofConstraints & constraints,
                     const dof_id_type their_dof_g =
                       parent_dof_indices[their_dof];
 
-                    const Real their_dof_value = FEInterface::shape(Dim-1,
-                                                                    side_fe_type,
+                    const Real their_dof_value = FEInterface::shape(side_fe_type,
                                                                     parent_side.get(),
                                                                     their_dof,
                                                                     mapped_point);
@@ -1058,6 +1073,11 @@ template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o)
 template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
 template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
 
+template <> unsigned int FE<0,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e, o); }
+template <> unsigned int FE<1,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e, o); }
+template <> unsigned int FE<2,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e, o); }
+template <> unsigned int FE<3,LAGRANGE>::n_dofs(const Elem * e, const Order o) { return lagrange_n_dofs(e, o); }
+
 
 // Do full-specialization for every dimension, instead
 // of explicit instantiation at the end of this function.
@@ -1066,6 +1086,11 @@ template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const ElemType t, const
 template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
 
+template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(e, o, n); }
+
 
 // Lagrange elements have no dofs per element
 // (just at the nodes)
@@ -1074,6 +1099,11 @@ template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const O
 template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 
+template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+
 // Lagrange FEMs are always C^0 continuous
 template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
 template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_lagrange_vec.C b/src/fe/fe_lagrange_vec.C
index b1d438bd7d0..8eb9b517c73 100644
--- a/src/fe/fe_lagrange_vec.C
+++ b/src/fe/fe_lagrange_vec.C
@@ -1309,17 +1309,26 @@ template <> RealGradient FE<3,L2_LAGRANGE_VEC>::shape_second_deriv(const Elem *
 
 // Do full-specialization for every dimension, instead
 // of explicit instantiation at the end of this function.
-// This could be macro-ified.
 template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,LAGRANGE>::n_dofs(t,o); }
 template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,LAGRANGE>::n_dofs(t,o); }
 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(t,o); }
 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(t,o); }
 
+template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,LAGRANGE>::n_dofs(e,o); }
+template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,LAGRANGE>::n_dofs(e,o); }
+template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,LAGRANGE>::n_dofs(e,o); }
+template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,LAGRANGE>::n_dofs(e,o); }
+
 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0,L2_LAGRANGE>::n_dofs(t,o); }
 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1,L2_LAGRANGE>::n_dofs(t,o); }
 template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 2*FE<2,L2_LAGRANGE>::n_dofs(t,o); }
 template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs(const ElemType t, const Order o) { return 3*FE<3,L2_LAGRANGE>::n_dofs(t,o); }
 
+template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0,L2_LAGRANGE>::n_dofs(e,o); }
+template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1,L2_LAGRANGE>::n_dofs(e,o); }
+template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 2*FE<2,L2_LAGRANGE>::n_dofs(e,o); }
+template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs(const Elem * e, const Order o) { return 3*FE<3,L2_LAGRANGE>::n_dofs(e,o); }
+
 
 // Do full-specialization for every dimension, instead
 // of explicit instantiation at the end of this function.
@@ -1328,11 +1337,21 @@ template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, co
 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 2*FE<2,LAGRANGE>::n_dofs_at_node(t,o,n); }
 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return 3*FE<3,LAGRANGE>::n_dofs_at_node(t,o,n); }
 
+template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,LAGRANGE>::n_dofs_at_node(e,o,n); }
+template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,LAGRANGE>::n_dofs_at_node(e,o,n); }
+template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 2*FE<2,LAGRANGE>::n_dofs_at_node(e,o,n); }
+template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return 3*FE<3,LAGRANGE>::n_dofs_at_node(e,o,n); }
+
 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
 
 // Lagrange elements have no dofs per element
 // (just at the nodes)
@@ -1341,12 +1360,22 @@ template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, con
 template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 
+template <> unsigned int FE<0,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<1,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<2,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<3,LAGRANGE_VEC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+
 // L2 Lagrange elements have all their dofs on the element
 template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,L2_LAGRANGE_VEC>::n_dofs(t, o); }
 template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,L2_LAGRANGE_VEC>::n_dofs(t, o); }
 template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<2,L2_LAGRANGE_VEC>::n_dofs(t, o); }
 template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<3,L2_LAGRANGE_VEC>::n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
+template <> unsigned int FE<1,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
+template <> unsigned int FE<2,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<2,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
+template <> unsigned int FE<3,L2_LAGRANGE_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<3,L2_LAGRANGE_VEC>::n_dofs(&e, o); }
+
 // Lagrange FEMs are always C^0 continuous
 template <> FEContinuity FE<0,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
 template <> FEContinuity FE<1,LAGRANGE_VEC>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_map.C b/src/fe/fe_map.C
index d66f5fbc92a..8f80633e5a5 100644
--- a/src/fe/fe_map.C
+++ b/src/fe/fe_map.C
@@ -2014,7 +2014,7 @@ Point FEMap::inverse_map (const unsigned int dim,
       // Make sure the point \p p on the reference element actually
       // is
 
-      if (!FEAbstract::on_reference_element(p, elem->type(), 2*tolerance))
+      if (!elem->on_reference_element(p, 2*tolerance))
         {
           libmesh_here();
           libMesh::err << "WARNING:  inverse_map of physical point "
diff --git a/src/fe/fe_monomial.C b/src/fe/fe_monomial.C
index 744674277cd..2ac3c119bb8 100644
--- a/src/fe/fe_monomial.C
+++ b/src/fe/fe_monomial.C
@@ -27,6 +27,13 @@
 namespace libMesh
 {
 
+unsigned int monomial_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return monomial_n_dofs(e->type(), o);
+}
+
+
 unsigned int monomial_n_dofs(const ElemType t, const Order o)
 {
   switch (o)
@@ -377,6 +384,11 @@ template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o)
 template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 
+template <> unsigned int FE<0,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+template <> unsigned int FE<1,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+template <> unsigned int FE<2,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+template <> unsigned int FE<3,MONOMIAL>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 // Monomials have no dofs at nodes, only element dofs.
 template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
@@ -384,12 +396,22 @@ template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Or
 template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 
+template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+
 
 // Full specialization of get_continuity() function for every dimension.
 template <> FEContinuity FE<0,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
diff --git a/src/fe/fe_monomial_vec.C b/src/fe/fe_monomial_vec.C
index bfabd9e97ca..ce444d114c5 100644
--- a/src/fe/fe_monomial_vec.C
+++ b/src/fe/fe_monomial_vec.C
@@ -694,80 +694,35 @@ FE<3, MONOMIAL_VEC>::shape_second_deriv(const Elem * elem,
 // Do full-specialization for every dimension, instead
 // of explicit instantiation at the end of this function.
 // This could be macro-ified.
-template <>
-unsigned int
-FE<0, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o)
-{
-  return FE<0, MONOMIAL>::n_dofs(t, o);
-}
-template <>
-unsigned int
-FE<1, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o)
-{
-  return FE<1, MONOMIAL>::n_dofs(t, o);
-}
-template <>
-unsigned int
-FE<2, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o)
-{
-  return 2 * FE<2, MONOMIAL>::n_dofs(t, o);
-}
-template <>
-unsigned int
-FE<3, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o)
-{
-  return 3 * FE<3, MONOMIAL>::n_dofs(t, o);
-}
-
-template <>
-unsigned int
-FE<0, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int)
-{
-  return 0;
-}
-template <>
-unsigned int
-FE<1, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int)
-{
-  return 0;
-}
-template <>
-unsigned int
-FE<2, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int)
-{
-  return 0;
-}
-template <>
-unsigned int
-FE<3, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int)
-{
-  return 0;
-}
-
-template <>
-unsigned int
-FE<0, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o)
-{
-  return FE<0, MONOMIAL>::n_dofs_per_elem(t, o);
-}
-template <>
-unsigned int
-FE<1, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o)
-{
-  return FE<1, MONOMIAL>::n_dofs_per_elem(t, o);
-}
-template <>
-unsigned int
-FE<2, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o)
-{
-  return 2 * FE<2, MONOMIAL>::n_dofs_per_elem(t, o);
-}
-template <>
-unsigned int
-FE<3, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o)
-{
-  return 3 * FE<3, MONOMIAL>::n_dofs_per_elem(t, o);
-}
+template <> unsigned int FE<0, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o) { return FE<0, MONOMIAL>::n_dofs(t, o); }
+template <> unsigned int FE<1, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o) { return FE<1, MONOMIAL>::n_dofs(t, o); }
+template <> unsigned int FE<2, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o) { return 2 * FE<2, MONOMIAL>::n_dofs(t, o); }
+template <> unsigned int FE<3, MONOMIAL_VEC>::n_dofs(const ElemType t, const Order o) { return 3 * FE<3, MONOMIAL>::n_dofs(t, o); }
+
+template <> unsigned int FE<0, MONOMIAL_VEC>::n_dofs(const Elem * e, const Order o) { return FE<0, MONOMIAL>::n_dofs(e, o); }
+template <> unsigned int FE<1, MONOMIAL_VEC>::n_dofs(const Elem * e, const Order o) { return FE<1, MONOMIAL>::n_dofs(e, o); }
+template <> unsigned int FE<2, MONOMIAL_VEC>::n_dofs(const Elem * e, const Order o) { return 2 * FE<2, MONOMIAL>::n_dofs(e, o); }
+template <> unsigned int FE<3, MONOMIAL_VEC>::n_dofs(const Elem * e, const Order o) { return 3 * FE<3, MONOMIAL>::n_dofs(e, o); }
+
+template <> unsigned int FE<0, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3, MONOMIAL_VEC>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
+
+template <> unsigned int FE<0, MONOMIAL_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1, MONOMIAL_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2, MONOMIAL_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3, MONOMIAL_VEC>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
+template <> unsigned int FE<0, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0, MONOMIAL>::n_dofs_per_elem(t, o); }
+template <> unsigned int FE<1, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1, MONOMIAL>::n_dofs_per_elem(t, o); }
+template <> unsigned int FE<2, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 2 * FE<2, MONOMIAL>::n_dofs_per_elem(t, o); }
+template <> unsigned int FE<3, MONOMIAL_VEC>::n_dofs_per_elem(const ElemType t, const Order o) { return 3 * FE<3, MONOMIAL>::n_dofs_per_elem(t, o); }
+
+template <> unsigned int FE<0, MONOMIAL_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0, MONOMIAL>::n_dofs_per_elem(e, o); }
+template <> unsigned int FE<1, MONOMIAL_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1, MONOMIAL>::n_dofs_per_elem(e, o); }
+template <> unsigned int FE<2, MONOMIAL_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 2 * FE<2, MONOMIAL>::n_dofs_per_elem(e, o); }
+template <> unsigned int FE<3, MONOMIAL_VEC>::n_dofs_per_elem(const Elem & e, const Order o) { return 3 * FE<3, MONOMIAL>::n_dofs_per_elem(e, o); }
 
 // Monomial FEMs are always C^0 continuous
 template <>
diff --git a/src/fe/fe_nedelec_one.C b/src/fe/fe_nedelec_one.C
index aea16bdcc2b..b29aff3234b 100644
--- a/src/fe/fe_nedelec_one.C
+++ b/src/fe/fe_nedelec_one.C
@@ -95,6 +95,7 @@ void nedelec_one_nodal_soln(const Elem * elem,
 } // nedelec_one_nodal_soln
 
 
+
 unsigned int nedelec_one_n_dofs(const ElemType t, const Order o)
 {
   libmesh_assert_greater (o, 0);
@@ -124,6 +125,15 @@ unsigned int nedelec_one_n_dofs(const ElemType t, const Order o)
 }
 
 
+
+unsigned int nedelec_one_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return nedelec_one_n_dofs(e->type(), o);
+}
+
+
+
 unsigned int nedelec_one_n_dofs_at_node(const ElemType t,
                                         const Order o,
                                         const unsigned int n)
@@ -254,6 +264,16 @@ unsigned int nedelec_one_n_dofs_at_node(const ElemType t,
 }
 
 
+
+unsigned int nedelec_one_n_dofs_at_node(const Elem & e,
+                                        const Order o,
+                                        const unsigned int n)
+{
+  return nedelec_one_n_dofs_at_node(e.type(), o, n);
+}
+
+
+
 unsigned int nedelec_one_n_dofs_per_elem(const ElemType t,
                                          const Order o)
 {
@@ -284,6 +304,15 @@ unsigned int nedelec_one_n_dofs_per_elem(const ElemType t,
 }
 
 
+
+unsigned int nedelec_one_n_dofs_per_elem(const Elem & e,
+                                         const Order o)
+{
+  return nedelec_one_n_dofs_per_elem(e.type(), o);
+}
+
+
+
 #ifdef LIBMESH_ENABLE_AMR
 void nedelec_one_compute_constraints (DofConstraints & /*constraints*/,
                                       DofMap & /*dof_map*/,
@@ -346,22 +375,36 @@ LIBMESH_FE_SIDE_NODAL_SOLN(NEDELEC_ONE)
 
 // Do full-specialization for every dimension, instead
 // of explicit instantiation at the end of this function.
-// This could be macro-ified.
 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
 
+template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const Elem *, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const Elem *, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const Elem * e, const Order o) { return nedelec_one_n_dofs(e, o); }
+template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const Elem * e, const Order o) { return nedelec_one_n_dofs(e, o); }
+
 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); }
 
+template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(e, o, n); }
+
 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType t, const Order o) { return nedelec_one_n_dofs_per_elem(t, o); }
 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType t, const Order o) { return nedelec_one_n_dofs_per_elem(t, o); }
 
+template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const Elem &, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const Elem &, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const Elem & e, const Order o) { return nedelec_one_n_dofs_per_elem(e, o); }
+template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const Elem & e, const Order o) { return nedelec_one_n_dofs_per_elem(e, o); }
+
 // Nedelec first type FEMs are always tangentially continuous
 template <> FEContinuity FE<0,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
 template <> FEContinuity FE<1,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
diff --git a/src/fe/fe_rational.C b/src/fe/fe_rational.C
index 0e8a162f263..2b698cabc6c 100644
--- a/src/fe/fe_rational.C
+++ b/src/fe/fe_rational.C
@@ -132,18 +132,33 @@ template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, cons
 template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return FE<2,_underlying_fe_family>::n_dofs(t, o); }
 template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return FE<3,_underlying_fe_family>::n_dofs(t, o); }
 
+template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<0,_underlying_fe_family>::n_dofs(e, o); }
+template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<1,_underlying_fe_family>::n_dofs(e, o); }
+template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<2,_underlying_fe_family>::n_dofs(e, o); }
+template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<3,_underlying_fe_family>::n_dofs(e, o); }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<2,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<3,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
 
+template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<2,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<3,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,_underlying_fe_family>::n_dofs_per_elem(t, o); }
 template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,_underlying_fe_family>::n_dofs_per_elem(t, o); }
 template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<2,_underlying_fe_family>::n_dofs_per_elem(t, o); }
 template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<3,_underlying_fe_family>::n_dofs_per_elem(t, o); }
 
+template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,_underlying_fe_family>::n_dofs_per_elem(e, o); }
+template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,_underlying_fe_family>::n_dofs_per_elem(e, o); }
+template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<2,_underlying_fe_family>::n_dofs_per_elem(e, o); }
+template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<3,_underlying_fe_family>::n_dofs_per_elem(e, o); }
+
 // Our current rational FEMs are C^0 continuous
 template <> FEContinuity FE<0,RATIONAL_BERNSTEIN>::get_continuity() const { return C_ZERO; }
 template <> FEContinuity FE<1,RATIONAL_BERNSTEIN>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_raviart.C b/src/fe/fe_raviart.C
index e8c677fe1f3..c02ee891a21 100644
--- a/src/fe/fe_raviart.C
+++ b/src/fe/fe_raviart.C
@@ -98,6 +98,7 @@ void raviart_thomas_nodal_soln(const Elem * elem,
 } // raviart_thomas_nodal_soln
 
 
+
 unsigned int raviart_thomas_n_dofs(const ElemType t, const Order o)
 {
   libmesh_assert_greater (o, 0);
@@ -121,6 +122,15 @@ unsigned int raviart_thomas_n_dofs(const ElemType t, const Order o)
 }
 
 
+
+unsigned int raviart_thomas_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return raviart_thomas_n_dofs(e->type(), o);
+}
+
+
+
 unsigned int raviart_thomas_n_dofs_at_node(const ElemType t,
                                            const Order o,
                                            const unsigned int n)
@@ -240,6 +250,16 @@ unsigned int raviart_thomas_n_dofs_at_node(const ElemType t,
 }
 
 
+
+unsigned int raviart_thomas_n_dofs_at_node(const Elem & e,
+                                           const Order o,
+                                           const unsigned int n)
+{
+  return raviart_thomas_n_dofs_at_node(e.type(), o, n);
+}
+
+
+
 unsigned int raviart_thomas_n_dofs_per_elem(const ElemType t,
                                             const Order o)
 {
@@ -264,6 +284,14 @@ unsigned int raviart_thomas_n_dofs_per_elem(const ElemType t,
 }
 
 
+
+unsigned int raviart_thomas_n_dofs_per_elem(const Elem & e,
+                                            const Order o)
+{
+  return raviart_thomas_n_dofs_per_elem(e.type(), o);
+}
+
+
 #ifdef LIBMESH_ENABLE_AMR
 void raviart_thomas_compute_constraints (DofConstraints & /*constraints*/,
                                          DofMap & /*dof_map*/,
@@ -364,31 +392,63 @@ template <> unsigned int FE<0,RAVIART_THOMAS>::n_dofs(const ElemType, const Orde
 template <> unsigned int FE<1,RAVIART_THOMAS>::n_dofs(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,RAVIART_THOMAS>::n_dofs(const ElemType t, const Order o) { return raviart_thomas_n_dofs(t, o); }
 template <> unsigned int FE<3,RAVIART_THOMAS>::n_dofs(const ElemType t, const Order o) { return raviart_thomas_n_dofs(t, o); }
+
+template <> unsigned int FE<0,RAVIART_THOMAS>::n_dofs(const Elem *, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,RAVIART_THOMAS>::n_dofs(const Elem *, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,RAVIART_THOMAS>::n_dofs(const Elem * e, const Order o) { return raviart_thomas_n_dofs(e, o); }
+template <> unsigned int FE<3,RAVIART_THOMAS>::n_dofs(const Elem * e, const Order o) { return raviart_thomas_n_dofs(e, o); }
+
 template <> unsigned int FE<0,L2_RAVIART_THOMAS>::n_dofs(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,L2_RAVIART_THOMAS>::n_dofs(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,L2_RAVIART_THOMAS>::n_dofs(const ElemType t, const Order o) { return raviart_thomas_n_dofs(t, o); }
 template <> unsigned int FE<3,L2_RAVIART_THOMAS>::n_dofs(const ElemType t, const Order o) { return raviart_thomas_n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_RAVIART_THOMAS>::n_dofs(const Elem *, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,L2_RAVIART_THOMAS>::n_dofs(const Elem *, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,L2_RAVIART_THOMAS>::n_dofs(const Elem * e, const Order o) { return raviart_thomas_n_dofs(e, o); }
+template <> unsigned int FE<3,L2_RAVIART_THOMAS>::n_dofs(const Elem * e, const Order o) { return raviart_thomas_n_dofs(e, o); }
+
 template <> unsigned int FE<0,RAVIART_THOMAS>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,RAVIART_THOMAS>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,RAVIART_THOMAS>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return raviart_thomas_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<3,RAVIART_THOMAS>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return raviart_thomas_n_dofs_at_node(t, o, n); }
+
+template <> unsigned int FE<0,RAVIART_THOMAS>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,RAVIART_THOMAS>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,RAVIART_THOMAS>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return raviart_thomas_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<3,RAVIART_THOMAS>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return raviart_thomas_n_dofs_at_node(e, o, n); }
+
 template <> unsigned int FE<0,L2_RAVIART_THOMAS>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,L2_RAVIART_THOMAS>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,L2_RAVIART_THOMAS>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,L2_RAVIART_THOMAS>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,L2_RAVIART_THOMAS>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,L2_RAVIART_THOMAS>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,L2_RAVIART_THOMAS>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,L2_RAVIART_THOMAS>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
 template <> unsigned int FE<0,RAVIART_THOMAS>::n_dofs_per_elem(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,RAVIART_THOMAS>::n_dofs_per_elem(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,RAVIART_THOMAS>::n_dofs_per_elem(const ElemType t, const Order o) { return raviart_thomas_n_dofs_per_elem(t, o); }
 template <> unsigned int FE<3,RAVIART_THOMAS>::n_dofs_per_elem(const ElemType t, const Order o) { return raviart_thomas_n_dofs_per_elem(t, o); }
 
+template <> unsigned int FE<0,RAVIART_THOMAS>::n_dofs_per_elem(const Elem &, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,RAVIART_THOMAS>::n_dofs_per_elem(const Elem &, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,RAVIART_THOMAS>::n_dofs_per_elem(const Elem & e, const Order o) { return raviart_thomas_n_dofs_per_elem(e, o); }
+template <> unsigned int FE<3,RAVIART_THOMAS>::n_dofs_per_elem(const Elem & e, const Order o) { return raviart_thomas_n_dofs_per_elem(e, o); }
+
 // L2 Raviart-Thomas elements have all their dofs per element
 template <> unsigned int FE<0,L2_RAVIART_THOMAS>::n_dofs_per_elem(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<1,L2_RAVIART_THOMAS>::n_dofs_per_elem(const ElemType, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> unsigned int FE<2,L2_RAVIART_THOMAS>::n_dofs_per_elem(const ElemType t, const Order o) { return n_dofs(t, o); }
 template <> unsigned int FE<3,L2_RAVIART_THOMAS>::n_dofs_per_elem(const ElemType t, const Order o) { return n_dofs(t, o); }
 
+template <> unsigned int FE<0,L2_RAVIART_THOMAS>::n_dofs_per_elem(const Elem &, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<1,L2_RAVIART_THOMAS>::n_dofs_per_elem(const Elem &, const Order) { RAVIART_LOW_D_ERROR_MESSAGE }
+template <> unsigned int FE<2,L2_RAVIART_THOMAS>::n_dofs_per_elem(const Elem & e, const Order o) { return n_dofs(&e, o); }
+template <> unsigned int FE<3,L2_RAVIART_THOMAS>::n_dofs_per_elem(const Elem & e, const Order o) { return n_dofs(&e, o); }
+
 // Raviart-Thomas FEMs are always normally continuous
 template <> FEContinuity FE<0,RAVIART_THOMAS>::get_continuity() const { RAVIART_LOW_D_ERROR_MESSAGE }
 template <> FEContinuity FE<1,RAVIART_THOMAS>::get_continuity() const { RAVIART_LOW_D_ERROR_MESSAGE }
diff --git a/src/fe/fe_scalar.C b/src/fe/fe_scalar.C
index 7ef1fcc5fc0..c23fab4eb8a 100644
--- a/src/fe/fe_scalar.C
+++ b/src/fe/fe_scalar.C
@@ -62,6 +62,11 @@ template <> unsigned int FE<1,SCALAR>::n_dofs(const ElemType, const Order o) { r
 template <> unsigned int FE<2,SCALAR>::n_dofs(const ElemType, const Order o) { return o; }
 template <> unsigned int FE<3,SCALAR>::n_dofs(const ElemType, const Order o) { return o; }
 
+template <> unsigned int FE<0,SCALAR>::n_dofs(const Elem *, const Order o) { return o; }
+template <> unsigned int FE<1,SCALAR>::n_dofs(const Elem *, const Order o) { return o; }
+template <> unsigned int FE<2,SCALAR>::n_dofs(const Elem *, const Order o) { return o; }
+template <> unsigned int FE<3,SCALAR>::n_dofs(const Elem *, const Order o) { return o; }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 // SCALARs have no dofs at nodes
 template <> unsigned int FE<0,SCALAR>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
@@ -69,6 +74,11 @@ template <> unsigned int FE<1,SCALAR>::n_dofs_at_node(const ElemType, const Orde
 template <> unsigned int FE<2,SCALAR>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,SCALAR>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,SCALAR>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,SCALAR>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,SCALAR>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,SCALAR>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 // SCALARs have no dofs per element
 template <> unsigned int FE<0,SCALAR>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
@@ -76,6 +86,11 @@ template <> unsigned int FE<1,SCALAR>::n_dofs_per_elem(const ElemType, const Ord
 template <> unsigned int FE<2,SCALAR>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 template <> unsigned int FE<3,SCALAR>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 
+template <> unsigned int FE<0,SCALAR>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<1,SCALAR>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<2,SCALAR>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<3,SCALAR>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+
 // Scalar FEMs are discontinuous
 template <> FEContinuity FE<0,SCALAR>::get_continuity() const { return DISCONTINUOUS; }
 template <> FEContinuity FE<1,SCALAR>::get_continuity() const { return DISCONTINUOUS; }
diff --git a/src/fe/fe_side_hierarchic.C b/src/fe/fe_side_hierarchic.C
index 2283465d7eb..87c4a7f5a6e 100644
--- a/src/fe/fe_side_hierarchic.C
+++ b/src/fe/fe_side_hierarchic.C
@@ -89,6 +89,7 @@ void side_hierarchic_side_nodal_soln
 }
 
 
+
 unsigned int side_hierarchic_n_dofs_at_node(const ElemType t,
                                             const Order o,
                                             const unsigned int n)
@@ -145,6 +146,15 @@ unsigned int side_hierarchic_n_dofs_at_node(const ElemType t,
 
 
 
+unsigned int side_hierarchic_n_dofs_at_node(const Elem & e,
+                                            const Order o,
+                                            const unsigned int n)
+{
+  return side_hierarchic_n_dofs_at_node(e.type(), o, n);
+}
+
+
+
 unsigned int side_hierarchic_n_dofs(const ElemType t, const Order o)
 {
   switch (t)
@@ -177,6 +187,13 @@ unsigned int side_hierarchic_n_dofs(const ElemType t, const Order o)
 } // side_hierarchic_n_dofs()
 
 
+
+unsigned int side_hierarchic_n_dofs(const Elem * e, const Order o)
+{
+  return side_hierarchic_n_dofs(e->type(), o);
+}
+
+
 } // anonymous namespace
 
 
@@ -247,18 +264,33 @@ template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const O
 template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
 template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return side_hierarchic_n_dofs(t, o); }
 
+template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
+template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
+template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
+template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs(const Elem * e, const Order o) { return side_hierarchic_n_dofs(e, o); }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(t, o, n); }
 
+template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return side_hierarchic_n_dofs_at_node(e, o, n); }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
 
+template <> unsigned int FE<0,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<1,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<2,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+template <> unsigned int FE<3,SIDE_HIERARCHIC>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
+
 // Side FEMs are discontinuous from side to side
 template <> FEContinuity FE<0,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
 template <> FEContinuity FE<1,SIDE_HIERARCHIC>::get_continuity() const { return SIDE_DISCONTINUOUS; }
diff --git a/src/fe/fe_subdivision_2D.C b/src/fe/fe_subdivision_2D.C
index b73aed51e78..56ec42cf5d3 100644
--- a/src/fe/fe_subdivision_2D.C
+++ b/src/fe/fe_subdivision_2D.C
@@ -954,12 +954,15 @@ void FE<2,SUBDIVISION>::inverse_map(const Elem *,
 
 // The number of dofs per subdivision element depends on the valence of its nodes, so it is non-static
 template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; }
+template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const Elem *, const Order) { libmesh_not_implemented(); return 0; }
 
 // Loop subdivision elements have only a single dof per node
 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; }
+template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 1; }
 
 // Subdivision FEMs have dofs only at the nodes
 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
+template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const Elem &, const Order) { return 0; }
 
 // Subdivision FEMs have dofs only at the nodes
 template <> void FE<2,SUBDIVISION>::dofs_on_side(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di, bool) { di.resize(0); }
diff --git a/src/fe/fe_szabab.C b/src/fe/fe_szabab.C
index 4073fd51a78..dccfbf4deff 100644
--- a/src/fe/fe_szabab.C
+++ b/src/fe/fe_szabab.C
@@ -109,7 +109,6 @@ void szabab_nodal_soln(const Elem * elem,
 
 
 
-
 unsigned int szabab_n_dofs(const ElemType t, const Order o)
 {
   switch (o)
@@ -326,6 +325,12 @@ unsigned int szabab_n_dofs(const ElemType t, const Order o)
 
 
 
+unsigned int szabab_n_dofs(const Elem * e, const Order o)
+{
+  libmesh_assert(e);
+  return szabab_n_dofs(e->type(), o);
+}
+
 
 
 unsigned int szabab_n_dofs_at_node(const ElemType t,
@@ -966,6 +971,15 @@ unsigned int szabab_n_dofs_at_node(const ElemType t,
 
 
 
+unsigned int szabab_n_dofs_at_node(const Elem & e,
+                                   const Order o,
+                                   const unsigned int n)
+{
+  return szabab_n_dofs_at_node(e.type(), o, n);
+}
+
+
+
 unsigned int szabab_n_dofs_per_elem(const ElemType t, const Order o)
 {
   switch (o)
@@ -1265,6 +1279,13 @@ unsigned int szabab_n_dofs_per_elem(const ElemType t, const Order o)
     }
 } // szabab_n_dofs_per_elem
 
+
+
+unsigned int szabab_n_dofs_per_elem(const Elem & e, const Order o)
+{
+  return szabab_n_dofs_per_elem(e.type(), o);
+}
+
 } // anonymous namespace
 
 
@@ -1279,18 +1300,33 @@ template <> unsigned int FE<1,SZABAB>::n_dofs(const ElemType t, const Order o) {
 template <> unsigned int FE<2,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
 template <> unsigned int FE<3,SZABAB>::n_dofs(const ElemType t, const Order o) { return szabab_n_dofs(t, o); }
 
+template <> unsigned int FE<0,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
+template <> unsigned int FE<1,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
+template <> unsigned int FE<2,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
+template <> unsigned int FE<3,SZABAB>::n_dofs(const Elem * e, const Order o) { return szabab_n_dofs(e, o); }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 template <> unsigned int FE<0,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<1,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<2,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
 template <> unsigned int FE<3,SZABAB>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(t, o, n); }
 
+template <> unsigned int FE<0,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<1,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<2,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
+template <> unsigned int FE<3,SZABAB>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return szabab_n_dofs_at_node(e, o, n); }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 template <> unsigned int FE<0,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
 template <> unsigned int FE<1,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
 template <> unsigned int FE<2,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
 template <> unsigned int FE<3,SZABAB>::n_dofs_per_elem(const ElemType t, const Order o) { return szabab_n_dofs_per_elem(t, o); }
 
+template <> unsigned int FE<0,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
+template <> unsigned int FE<1,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
+template <> unsigned int FE<2,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
+template <> unsigned int FE<3,SZABAB>::n_dofs_per_elem(const Elem & e, const Order o) { return szabab_n_dofs_per_elem(e, o); }
+
 // Szabab FEMs are C^0 continuous
 template <> FEContinuity FE<0,SZABAB>::get_continuity() const { return C_ZERO; }
 template <> FEContinuity FE<1,SZABAB>::get_continuity() const { return C_ZERO; }
diff --git a/src/fe/fe_xyz.C b/src/fe/fe_xyz.C
index aff571cf8f9..367cd8c028b 100644
--- a/src/fe/fe_xyz.C
+++ b/src/fe/fe_xyz.C
@@ -102,7 +102,7 @@ void xyz_nodal_soln(const Elem * elem,
 
 template <unsigned int Dim>
 void FEXYZ<Dim>::init_shape_functions(const std::vector<Point> & qp,
-                                      const Elem * libmesh_dbg_var(elem))
+                                      const Elem * elem)
 {
   libmesh_assert(elem);
 
@@ -119,8 +119,8 @@ void FEXYZ<Dim>::init_shape_functions(const std::vector<Point> & qp,
   // Number of shape functions in the finite element approximation
   // space.
   const unsigned int n_approx_shape_functions =
-    this->n_shape_functions(this->get_type(),
-                            this->get_order());
+    this->n_dofs(elem,
+                 this->get_order());
 
   // resize the vectors to hold current data
   // Phi are the shape functions used for the FE approximation
@@ -380,6 +380,11 @@ template <> unsigned int FE<1,XYZ>::n_dofs(const ElemType t, const Order o) { re
 template <> unsigned int FE<2,XYZ>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<3,XYZ>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 
+template <> unsigned int FE<0,XYZ>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+template <> unsigned int FE<1,XYZ>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+template <> unsigned int FE<2,XYZ>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+template <> unsigned int FE<3,XYZ>::n_dofs(const Elem * e, const Order o) { return monomial_n_dofs(e, o); }
+
 // Full specialization of n_dofs_at_node() function for every dimension.
 // XYZ FEMs have no dofs at nodes, only element dofs.
 template <> unsigned int FE<0,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
@@ -387,12 +392,22 @@ template <> unsigned int FE<1,XYZ>::n_dofs_at_node(const ElemType, const Order,
 template <> unsigned int FE<2,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 template <> unsigned int FE<3,XYZ>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
 
+template <> unsigned int FE<0,XYZ>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<1,XYZ>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<2,XYZ>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+template <> unsigned int FE<3,XYZ>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { return 0; }
+
 // Full specialization of n_dofs_per_elem() function for every dimension.
 template <> unsigned int FE<0,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<1,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<2,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 template <> unsigned int FE<3,XYZ>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
 
+template <> unsigned int FE<0,XYZ>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+template <> unsigned int FE<1,XYZ>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+template <> unsigned int FE<2,XYZ>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+template <> unsigned int FE<3,XYZ>::n_dofs_per_elem(const Elem & e, const Order o) { return monomial_n_dofs(&e, o); }
+
 // Full specialization of get_continuity() function for every dimension.
 template <> FEContinuity FE<0,XYZ>::get_continuity() const { return DISCONTINUOUS; }
 template <> FEContinuity FE<1,XYZ>::get_continuity() const { return DISCONTINUOUS; }
diff --git a/src/fe/fe_xyz_boundary.C b/src/fe/fe_xyz_boundary.C
index 1d76f33d6e2..8ec6e460237 100644
--- a/src/fe/fe_xyz_boundary.C
+++ b/src/fe/fe_xyz_boundary.C
@@ -146,8 +146,7 @@ void FEXYZ<Dim>::compute_face_values(const Elem * elem,
   // Number of shape functions in the finite element approximation
   // space.
   const unsigned int n_approx_shape_functions =
-    this->n_shape_functions(this->get_type(),
-                            this->get_order());
+    this->n_dofs(elem, this->get_order());
 
   // Resize the shape functions and their gradients
   this->phi.resize    (n_approx_shape_functions);
diff --git a/src/fe/inf_fe.C b/src/fe/inf_fe.C
index dacf84411ad..c85915a35d0 100644
--- a/src/fe/inf_fe.C
+++ b/src/fe/inf_fe.C
@@ -510,7 +510,8 @@ void InfFE<Dim,T_radial,T_map>::init_shape_functions(const std::vector<Point> &
   // initialize most of the things related to physical approximation
   unsigned int n_base_approx_shape_functions;
   if (Dim > 1)
-    n_base_approx_shape_functions = base_fe->n_shape_functions();
+    n_base_approx_shape_functions =
+      FEInterface::n_dofs(base_fe->get_fe_type(), base_elem.get());
   else
     n_base_approx_shape_functions = 1;
 
@@ -819,7 +820,7 @@ void InfFE<Dim,T_radial,T_map>::compute_shape_functions(const Elem * inf_elem,
 
 
   _n_total_approx_sf = InfFERadial::n_dofs(fe_type.radial_order) *
-    base_fe->n_shape_functions();
+    FEInterface::n_dofs(base_fe->get_fe_type(), base_elem.get());
 
 
 
diff --git a/src/fe/inf_fe_static.C b/src/fe/inf_fe_static.C
index bfa7b42ada7..73463022c80 100644
--- a/src/fe/inf_fe_static.C
+++ b/src/fe/inf_fe_static.C
@@ -62,6 +62,7 @@ bool InfFE<Dim,T_radial,T_map>::_warned_for_dshape     = false;
 
 // ------------------------------------------------------------
 // InfFE static class members
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
 unsigned int InfFE<Dim,T_radial,T_map>::n_dofs (const FEType & fet,
                                                 const ElemType inf_elem_type)
@@ -76,6 +77,7 @@ unsigned int InfFE<Dim,T_radial,T_map>::n_dofs (const FEType & fet,
   else
     return InfFERadial::n_dofs(fet.radial_order);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -97,13 +99,13 @@ unsigned int InfFE<Dim,T_radial,T_map>::n_dofs(const FEType & fet,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
 unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_at_node (const FEType & fet,
                                                         const ElemType inf_elem_type,
                                                         const unsigned int n)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
 
@@ -123,6 +125,7 @@ unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_at_node (const FEType & fet,
   else
     return InfFERadial::n_dofs_at_node(fet.radial_order, n_radial);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -148,12 +151,12 @@ unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_at_node (const FEType & fet,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
 unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_per_elem (const FEType & fet,
                                                          const ElemType inf_elem_type)
 {
-  // TODO:
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   const ElemType base_et (InfFEBase::get_elem_type(inf_elem_type));
 
@@ -163,6 +166,7 @@ unsigned int InfFE<Dim,T_radial,T_map>::n_dofs_per_elem (const FEType & fet,
   else
     return InfFERadial::n_dofs_per_elem(fet.radial_order);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -219,14 +223,14 @@ void InfFE<Dim,T_radial,T_map>::nodal_soln(const FEType & /* fet */,
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
 Real InfFE<Dim,T_radial,T_map>::shape(const FEType & fet,
                                       const ElemType inf_elem_type,
                                       const unsigned int i,
                                       const Point & p)
 {
-  // TODO - if possible, not sure if we can easily fully remove this function.
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   libmesh_assert_not_equal_to (Dim, 0);
 
@@ -260,9 +264,7 @@ Real InfFE<Dim,T_radial,T_map>::shape(const FEType & fet,
     return InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
       * InfFERadial::decay(Dim,v);
 }
-
-
-
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -322,6 +324,7 @@ Real InfFE<Dim,T_radial,T_map>::shape(const FEType fet,
 }
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
 Real InfFE<Dim,T_radial,T_map>::shape_deriv (const FEType & fe_t,
                                              const ElemType inf_elem_type,
@@ -329,8 +332,7 @@ Real InfFE<Dim,T_radial,T_map>::shape_deriv (const FEType & fe_t,
                                              const unsigned int j,
                                              const Point & p)
 {
-  // TODO - if possible, not sure if we can easily fully remove this function.
-  // libmesh_deprecated();
+  libmesh_deprecated();
 
   libmesh_assert_not_equal_to (Dim, 0);
   libmesh_assert_greater (Dim,j);
@@ -367,6 +369,7 @@ Real InfFE<Dim,T_radial,T_map>::shape_deriv (const FEType & fe_t,
     * InfFE<Dim,T_radial,T_map>::eval(v, o_radial, i_radial)
     * InfFERadial::decay(Dim,v);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
@@ -948,27 +951,6 @@ void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
                                                        unsigned int & base_shape,
                                                        unsigned int & radial_shape)
 {
-  // (Temporarily) call version of this function taking an
-  // ElemType. Eventually there should only be one version of this
-  // function that takes an Elem*.
-  compute_shape_indices(fet, inf_elem->type(), i, base_shape, radial_shape);
-}
-
-
-
-template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
-void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
-                                                       const ElemType inf_elem_type,
-                                                       const unsigned int i,
-                                                       unsigned int & base_shape,
-                                                       unsigned int & radial_shape)
-{
-  // TODO: eventually figure out a way to deprecated this
-  // function. Note that we can't go the other way around and have
-  // this function call the Elem* version because there's not really a
-  // clean way to create the required Elem object on the fly...
-  // libmesh_deprecated();
-
   // An example is provided:  the numbers in comments refer to
   // a fictitious InfHex18.  The numbers are chosen as exemplary
   // values.  There is currently no base approximation that
@@ -995,6 +977,220 @@ void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
   const unsigned int radial_order       = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
   const unsigned int radial_order_p_one = radial_order+1;                                          // 5
 
+  std::unique_ptr<const Elem> base_elem = InfFEBase::build_elem(inf_elem);                         // QUAD9
+
+  // assume that the number of dof is the same for all vertices
+  unsigned int n_base_vertices         = libMesh::invalid_uint;                                    // 4
+  const unsigned int n_base_vertex_dof = FEInterface::n_dofs_at_node  (fet, base_elem.get(), 0);   // 2
+
+  unsigned int n_base_side_nodes       = libMesh::invalid_uint;                                    // 4
+  unsigned int n_base_side_dof         = libMesh::invalid_uint;                                    // 3
+
+  unsigned int n_base_face_nodes       = libMesh::invalid_uint;                                    // 1
+  unsigned int n_base_face_dof         = libMesh::invalid_uint;                                    // 5
+
+  const unsigned int n_base_elem_dof   = FEInterface::n_dofs_per_elem (fet, base_elem.get());      // 9
+
+  const ElemType inf_elem_type = inf_elem->type();
+
+  switch (inf_elem_type)
+    {
+    case INFEDGE2:
+      {
+        n_base_vertices   = 1;
+        n_base_side_nodes = 0;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = 0;
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    case INFQUAD4:
+      {
+        n_base_vertices   = 2;
+        n_base_side_nodes = 0;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = 0;
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    case INFQUAD6:
+      {
+        n_base_vertices   = 2;
+        n_base_side_nodes = 1;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    case INFHEX8:
+      {
+        n_base_vertices   = 4;
+        n_base_side_nodes = 0;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = 0;
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    case INFHEX16:
+      {
+        n_base_vertices   = 4;
+        n_base_side_nodes = 4;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    case INFHEX18:
+      {
+        n_base_vertices   = 4;
+        n_base_side_nodes = 4;
+        n_base_face_nodes = 1;
+        n_base_side_dof   = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
+        n_base_face_dof   = FEInterface::n_dofs_at_node (fet, base_elem.get(), 8);
+        break;
+      }
+
+
+    case INFPRISM6:
+      {
+        n_base_vertices   = 3;
+        n_base_side_nodes = 0;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = 0;
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    case INFPRISM12:
+      {
+        n_base_vertices   = 3;
+        n_base_side_nodes = 3;
+        n_base_face_nodes = 0;
+        n_base_side_dof   = FEInterface::n_dofs_at_node (fet, base_elem.get(), n_base_vertices);
+        n_base_face_dof   = 0;
+        break;
+      }
+
+    default:
+      libmesh_error_msg("Unrecognized inf_elem type = " << Utility::enum_to_string(inf_elem_type));
+    }
+
+
+  {
+    // these are the limits describing the intervals where the shape function lies
+    const unsigned int n_dof_at_base_vertices = n_base_vertices*n_base_vertex_dof;                 // 8
+    const unsigned int n_dof_at_all_vertices  = n_dof_at_base_vertices*radial_order_p_one;         // 40
+
+    const unsigned int n_dof_at_base_sides    = n_base_side_nodes*n_base_side_dof;                 // 12
+    const unsigned int n_dof_at_all_sides     = n_dof_at_base_sides*radial_order_p_one;            // 60
+
+    const unsigned int n_dof_at_base_face     = n_base_face_nodes*n_base_face_dof;                 // 5
+    const unsigned int n_dof_at_all_faces     = n_dof_at_base_face*radial_order_p_one;             // 25
+
+
+    // start locating the shape function
+    if (i < n_dof_at_base_vertices)                                              // range of i: 0..7
+      {
+        // belongs to vertex in the base
+        radial_shape = 0;
+        base_shape   = i;
+      }
+
+    else if (i < n_dof_at_all_vertices)                                          // range of i: 8..39
+      {
+        /* belongs to vertex in the outer shell
+         *
+         * subtract the number of dof already counted,
+         * so that i_offset contains only the offset for the base
+         */
+        const unsigned int i_offset = i - n_dof_at_base_vertices;                // 0..31
+
+        // first the radial dof are counted, then the base dof
+        radial_shape = (i_offset % radial_order) + 1;
+        base_shape   = i_offset / radial_order;
+      }
+
+    else if (i < n_dof_at_all_vertices+n_dof_at_base_sides)                      // range of i: 40..51
+      {
+        // belongs to base, is a side node
+        radial_shape = 0;
+        base_shape = i - radial_order * n_dof_at_base_vertices;                  //  8..19
+      }
+
+    else if (i < n_dof_at_all_vertices+n_dof_at_all_sides)                       // range of i: 52..99
+      {
+        // belongs to side node in the outer shell
+        const unsigned int i_offset = i - (n_dof_at_all_vertices
+                                           + n_dof_at_base_sides);               // 0..47
+        radial_shape = (i_offset % radial_order) + 1;
+        base_shape   = (i_offset / radial_order) + n_dof_at_base_vertices;
+      }
+
+    else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_base_face)    // range of i: 100..104
+      {
+        // belongs to the node in the base face
+        radial_shape = 0;
+        base_shape = i - radial_order*(n_dof_at_base_vertices
+                                       + n_dof_at_base_sides);                   //  20..24
+      }
+
+    else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces)    // range of i: 105..124
+      {
+        // belongs to the node in the outer face
+        const unsigned int i_offset = i - (n_dof_at_all_vertices
+                                           + n_dof_at_all_sides
+                                           + n_dof_at_base_face);                // 0..19
+        radial_shape = (i_offset % radial_order) + 1;
+        base_shape   = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides;
+      }
+
+    else if (i < n_dof_at_all_vertices+n_dof_at_all_sides+n_dof_at_all_faces+n_base_elem_dof)      // range of i: 125..133
+      {
+        // belongs to the base and is an element associated shape
+        radial_shape = 0;
+        base_shape = i - (n_dof_at_all_vertices
+                          + n_dof_at_all_sides
+                          + n_dof_at_all_faces);                                 // 0..8
+      }
+
+    else                                                                         // range of i: 134..169
+      {
+        libmesh_assert_less (i, n_dofs(fet, inf_elem));
+        // belongs to the outer shell and is an element associated shape
+        const unsigned int i_offset = i - (n_dof_at_all_vertices
+                                           + n_dof_at_all_sides
+                                           + n_dof_at_all_faces
+                                           + n_base_elem_dof);                   // 0..19
+        radial_shape = (i_offset % radial_order) + 1;
+        base_shape   = (i_offset / radial_order) + n_dof_at_base_vertices + n_dof_at_base_sides + n_dof_at_base_face;
+      }
+  }
+
+  return;
+}
+
+
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
+template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
+void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
+                                                       const ElemType inf_elem_type,
+                                                       const unsigned int i,
+                                                       unsigned int & base_shape,
+                                                       unsigned int & radial_shape)
+{
+  // This is basically cut-and-paste copy of the non-deprecated code
+  // above, and should be removed eventually.
+  libmesh_deprecated();
+
+  const unsigned int radial_order       = static_cast<unsigned int>(fet.radial_order.get_order()); // 4
+  const unsigned int radial_order_p_one = radial_order+1;                                          // 5
+
   const ElemType base_elem_type           (InfFEBase::get_elem_type(inf_elem_type));               // QUAD9
 
   // assume that the number of dof is the same for all vertices
@@ -1190,6 +1386,7 @@ void InfFE<Dim,T_radial,T_map>::compute_shape_indices (const FEType & fet,
 
   return;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 #ifdef LIBMESH_ENABLE_AMR
@@ -1454,8 +1651,7 @@ void InfFE<Dim, T_radial, T_map>::inf_compute_constraints (DofConstraints & cons
                        parent_base_dof != n_base_dofs; parent_base_dof++)
                     {
 
-                      const Real parent_base_dof_value = FEInterface::shape(Dim-1,
-                                                                            fe_type,
+                      const Real parent_base_dof_value = FEInterface::shape(fe_type,
                                                                             parent_base.get(),
                                                                             parent_base_dof,
                                                                             mapped_point);
@@ -1519,48 +1715,51 @@ void InfFE<Dim, T_radial, T_map>::inf_compute_constraints (DofConstraints & cons
 //#include "libmesh/inf_fe_instantiate_2D.h"
 //#include "libmesh/inf_fe_instantiate_3D.h"
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs(const FEType &, const ElemType));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs(const FEType &, const ElemType));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs(const FEType &, const ElemType));
-INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
-INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
-INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const ElemType));
-INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
-INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
-INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const ElemType, const unsigned int));
-INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
-INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
-INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const ElemType, const unsigned int, unsigned int &, unsigned int &));
-INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
-INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
-INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_node_indices(const ElemType, const unsigned int, unsigned int &, unsigned int &));
+INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
+INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
+INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
+INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
+INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
+INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
+#endif // LIBMESH_ENABLE_DEPRECATED
+
+INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
+INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
+INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs(const FEType &, const Elem*));
+INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
+INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
+INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_per_elem(const FEType &, const Elem *));
+INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
+INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
+INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, unsigned int, n_dofs_at_node(const FEType &, const Elem *, const unsigned int));
+INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
+INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
+INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, compute_shape_indices(const FEType &, const Elem *, const unsigned int, unsigned int &, unsigned int &));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const Elem *, const unsigned int, const Point & p));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType, const Elem *, const unsigned int, const Point &, const bool));
-INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
-INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
-INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape(const FEType &, const ElemType, const unsigned int, const Point &));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType &, const Elem *, const unsigned int, const unsigned int, const Point &));
-INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
-INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
-INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType &, const ElemType, const unsigned int, const unsigned int, const Point &));
 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Real, shape_deriv(const FEType, const Elem *, const unsigned int, const unsigned int, const Point &, const bool));
diff --git a/src/geom/cell_hex.C b/src/geom/cell_hex.C
index 70966cb5865..27958e7e7c0 100644
--- a/src/geom/cell_hex.C
+++ b/src/geom/cell_hex.C
@@ -574,6 +574,24 @@ std::pair<Real, Real> Hex::qual_bounds (const ElemQuality q) const
 
 
 
+bool Hex::on_reference_element(const Point & p,
+                               const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  const Real & zeta = p(2);
+
+  // The reference hexahedral element is [-1,1]^3.
+  return ((xi   >= -1.-eps) &&
+          (xi   <=  1.+eps) &&
+          (eta  >= -1.-eps) &&
+          (eta  <=  1.+eps) &&
+          (zeta >= -1.-eps) &&
+          (zeta <=  1.+eps));
+}
+
+
+
 const unsigned short int Hex::_second_order_vertex_child_number[27] =
   {
     99,99,99,99,99,99,99,99, // Vertices
diff --git a/src/geom/cell_inf_hex.C b/src/geom/cell_inf_hex.C
index f6461eae6bf..02e62cb17d8 100644
--- a/src/geom/cell_inf_hex.C
+++ b/src/geom/cell_inf_hex.C
@@ -533,9 +533,29 @@ bool InfHex::contains_point (const Point & p, Real tol) const
   const Point mapped_point = InfFEMap::inverse_map(dim(), this, p,
                                                    tol, false);
 
-  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
+  return this->on_reference_element(mapped_point, tol);
 }
 
+
+
+bool InfHex::on_reference_element(const Point & p,
+                                  const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  const Real & zeta = p(2);
+
+  // The reference infhex is a [-1,1]^3.
+  return ((xi   >= -1.-eps) &&
+          (xi   <=  1.+eps) &&
+          (eta  >= -1.-eps) &&
+          (eta  <=  1.+eps) &&
+          (zeta >= -1.-eps) &&
+          (zeta <=  1.+eps));
+}
+
+
+
 } // namespace libMesh
 
 
diff --git a/src/geom/cell_inf_prism.C b/src/geom/cell_inf_prism.C
index dfa1538b544..a7e059e0e04 100644
--- a/src/geom/cell_inf_prism.C
+++ b/src/geom/cell_inf_prism.C
@@ -329,7 +329,7 @@ bool InfPrism::contains_point (const Point & p, Real tol) const
   const Point mapped_point = InfFEMap::inverse_map(dim(), this, p,
                                                    tol, false);
 
-  return FEInterface::on_reference_element(mapped_point, this->type(), tol);
+  return this->on_reference_element(mapped_point, tol);
 }
 
 std::vector<unsigned>
@@ -371,6 +371,25 @@ InfPrism::edges_adjacent_to_node(const unsigned int n) const
 }
 
 
+
+bool InfPrism::on_reference_element(const Point & p,
+                                    const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  const Real & zeta = p(2);
+
+  // inside the reference triangle with zeta in [-1,1]
+  return ((xi   >=  0.-eps) &&
+          (eta  >=  0.-eps) &&
+          (zeta >= -1.-eps) &&
+          (zeta <=  1.+eps) &&
+          ((xi + eta) <= 1.+eps));
+}
+
+
+
+
 } // namespace libMesh
 
 
diff --git a/src/geom/cell_prism.C b/src/geom/cell_prism.C
index 00d0b42ee9d..6e48e77af9f 100644
--- a/src/geom/cell_prism.C
+++ b/src/geom/cell_prism.C
@@ -303,6 +303,23 @@ Prism::edges_adjacent_to_node(const unsigned int n) const
 
 
 
+bool Prism::on_reference_element(const Point & p,
+                                 const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  const Real & zeta = p(2);
+
+  // inside the reference triangle with zeta in [-1,1]
+  return ((xi   >=  0.-eps) &&
+          (eta  >=  0.-eps) &&
+          (zeta >= -1.-eps) &&
+          (zeta <=  1.+eps) &&
+          ((xi + eta) <= 1.+eps));
+}
+
+
+
 const unsigned short int Prism::_second_order_vertex_child_number[18] =
   {
     99,99,99,99,99,99, // Vertices
diff --git a/src/geom/cell_pyramid.C b/src/geom/cell_pyramid.C
index 61e8f02ba78..124e782613e 100644
--- a/src/geom/cell_pyramid.C
+++ b/src/geom/cell_pyramid.C
@@ -297,4 +297,27 @@ unsigned int Pyramid::local_singular_node(const Point & p, const Real tol) const
 }
 
 
+bool Pyramid::on_reference_element(const Point & p,
+                                   const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  const Real & zeta = p(2);
+
+  // Check that the point is on the same side of all the faces
+  // by testing whether:
+  //
+  // n_i.(x - x_i) <= 0
+  //
+  // for each i, where:
+  //   n_i is the outward normal of face i,
+  //   x_i is a point on face i.
+  return ((-eta - 1. + zeta <= 0.+eps) &&
+          (  xi - 1. + zeta <= 0.+eps) &&
+          ( eta - 1. + zeta <= 0.+eps) &&
+          ( -xi - 1. + zeta <= 0.+eps) &&
+          (            zeta >= 0.-eps));
+}
+
+
 } // namespace libMesh
diff --git a/src/geom/cell_tet.C b/src/geom/cell_tet.C
index 793176298f4..02bf693ab4e 100644
--- a/src/geom/cell_tet.C
+++ b/src/geom/cell_tet.C
@@ -324,4 +324,22 @@ std::pair<Real, Real> Tet::qual_bounds (const ElemQuality q) const
   return bounds;
 }
 
+
+
+bool Tet::on_reference_element(const Point & p,
+                               const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  const Real & zeta = p(2);
+        // The reference tetrahedral is isosceles
+        // and is bound by xi=0, eta=0, zeta=0,
+        // and xi+eta+zeta=1.
+  return ((xi   >= 0.-eps) &&
+          (eta  >= 0.-eps) &&
+          (zeta >= 0.-eps) &&
+          ((xi + eta + zeta) <= 1.+eps));
+}
+
+
 } // namespace libMesh
diff --git a/src/geom/edge.C b/src/geom/edge.C
index bb383cd7abc..39fc0c214b6 100644
--- a/src/geom/edge.C
+++ b/src/geom/edge.C
@@ -175,4 +175,11 @@ Edge::is_flipped() const
 }
 
 
+
+bool Edge::on_reference_element(const Point & p,
+                                const Real eps) const
+{
+  return ((p(0) >= -1-eps) && (p(0) <= 1+eps));
+}
+
 } // namespace libMesh
diff --git a/src/geom/elem.C b/src/geom/elem.C
index 932679e5695..80792f0413a 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -400,6 +400,7 @@ const Elem * Elem::reference_elem () const
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 Point Elem::centroid() const
 {
   libmesh_do_once(libMesh::err
@@ -412,6 +413,9 @@ Point Elem::centroid() const
 
   return Elem::vertex_average();
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
+
+
 
 Point Elem::true_centroid() const
 {
@@ -739,12 +743,14 @@ unsigned int Elem::which_side_am_i (const Elem * e) const
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 unsigned int Elem::which_node_am_i(unsigned int side,
                                    unsigned int side_node) const
 {
   libmesh_deprecated();
   return local_side_node(side, side_node);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
@@ -2684,14 +2690,12 @@ bool Elem::point_test(const Point & p, Real box_tol, Real map_tol) const
       // If dist is larger than some fraction of the tolerance, then return false.
       // This can happen when e.g. a 2D element is living in 3D, and
       // FEMap::inverse_map() maps p onto the projection of the element,
-      // effectively "tricking" FEInterface::on_reference_element().
+      // effectively "tricking" on_reference_element().
       if (dist > this->hmax() * map_tol)
         return false;
     }
 
-
-
-  return FEInterface::on_reference_element(mapped_point, this->type(), map_tol);
+  return this->on_reference_element(mapped_point, map_tol);
 }
 
 
diff --git a/src/geom/face_inf_quad.C b/src/geom/face_inf_quad.C
index 6e6db75584c..572ceaf0cf0 100644
--- a/src/geom/face_inf_quad.C
+++ b/src/geom/face_inf_quad.C
@@ -305,6 +305,23 @@ std::pair<Real, Real> InfQuad::qual_bounds (const ElemQuality q) const
   return bounds;
 }
 
+
+
+bool InfQuad::on_reference_element(const Point & p,
+                                   const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+
+  // The reference infquad is [-1,1]^2.
+  return ((xi  >= -1.-eps) &&
+          (xi  <=  1.+eps) &&
+          (eta >= -1.-eps) &&
+          (eta <=  1.+eps));
+}
+
+
+
 } // namespace libMesh
 
 
diff --git a/src/geom/face_inf_quad4.C b/src/geom/face_inf_quad4.C
index a16418c8a44..2c7a7195f55 100644
--- a/src/geom/face_inf_quad4.C
+++ b/src/geom/face_inf_quad4.C
@@ -148,7 +148,7 @@ bool InfQuad4::contains_point (const Point & p, Real tol) const
       const Point mapped_point = InfFEMap::inverse_map(dim(), this, p,
                                                        tol, false);
 
-      return FEInterface::on_reference_element(mapped_point, this->type(), tol);
+      return this->on_reference_element(mapped_point, tol);
     }
 }
 
diff --git a/src/geom/face_quad.C b/src/geom/face_quad.C
index c2ef64b67a7..cf6125d11a1 100644
--- a/src/geom/face_quad.C
+++ b/src/geom/face_quad.C
@@ -562,6 +562,20 @@ std::pair<Real, Real> Quad::qual_bounds (const ElemQuality q) const
 
 
 
+bool Quad::on_reference_element(const Point & p,
+                                const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+
+  // The reference quadrilateral element is [-1,1]^2.
+  return ((xi  >= -1.-eps) &&
+          (xi  <=  1.+eps) &&
+          (eta >= -1.-eps) &&
+          (eta <=  1.+eps));
+}
+
+
 
 const unsigned short int Quad::_second_order_adjacent_vertices[4][2] =
   {
diff --git a/src/geom/face_tri.C b/src/geom/face_tri.C
index b2aaee42973..34418e8e996 100644
--- a/src/geom/face_tri.C
+++ b/src/geom/face_tri.C
@@ -371,4 +371,16 @@ std::pair<Real, Real> Tri::qual_bounds (const ElemQuality q) const
   return bounds;
 }
 
+
+bool Tri::on_reference_element(const Point & p,
+                               const Real eps) const
+{
+  const Real & xi = p(0);
+  const Real & eta = p(1);
+  return ((xi  >= 0.-eps) &&
+          (eta >= 0.-eps) &&
+          ((xi + eta) <= 1.+eps));
+}
+
+
 } // namespace libMesh
diff --git a/src/geom/node_elem.C b/src/geom/node_elem.C
index 48e316cb6dc..38ce47624dc 100644
--- a/src/geom/node_elem.C
+++ b/src/geom/node_elem.C
@@ -65,6 +65,23 @@ bool NodeElem::close_to_point(const Point & p, Real tol) const
   return this->contains_point(p, tol);
 }
 
+
+bool NodeElem::on_reference_element(const Point & p,
+                                    const Real eps) const
+{
+  // Debatably `return true;` would make sense here; we shouldn't ever
+  // pass p!=0 in.
+  return (std::abs(p(0)) < eps
+#if LIBMESH_DIM > 1
+          && std::abs(p(1)) < eps
+#endif
+#if LIBMESH_DIM > 2
+          && std::abs(p(2)) < eps
+#endif
+         );
+}
+
+
 #ifdef LIBMESH_ENABLE_AMR
 
 const Real NodeElem::_embedding_matrix[1][1][1] =
diff --git a/src/mesh/dyna_io.C b/src/mesh/dyna_io.C
index 83c4f2d5ba1..c566239e2ab 100644
--- a/src/mesh/dyna_io.C
+++ b/src/mesh/dyna_io.C
@@ -740,12 +740,14 @@ void DynaIO::add_spline_constraints(DofMap &,
 }
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void DynaIO::clear_spline_nodes()
 {
   libmesh_deprecated();
 
   MeshTools::clear_spline_nodes(MeshInput<MeshBase>::mesh());
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/mesh/mesh_base.C b/src/mesh/mesh_base.C
index c7bfded07e4..a193f70fcda 100644
--- a/src/mesh/mesh_base.C
+++ b/src/mesh/mesh_base.C
@@ -737,6 +737,7 @@ void MeshBase::remove_orphaned_nodes ()
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
 {
   libmesh_deprecated();
@@ -768,6 +769,9 @@ void MeshBase::prepare_for_use (const bool skip_renumber_nodes_and_elements)
 
   this->prepare_for_use();
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
+
+
 
 void MeshBase::prepare_for_use ()
 {
@@ -1699,12 +1703,14 @@ subdomain_id_type MeshBase::get_id_by_name(std::string_view name) const
   return Elem::invalid_subdomain_id;
 }
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void MeshBase::cache_elem_dims()
 {
   libmesh_deprecated();
 
   this->cache_elem_data();
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 void MeshBase::cache_elem_data()
 {
diff --git a/src/mesh/mesh_communication.C b/src/mesh/mesh_communication.C
index bf15a9c364b..39a74779af5 100644
--- a/src/mesh/mesh_communication.C
+++ b/src/mesh/mesh_communication.C
@@ -349,6 +349,7 @@ void connect_children(const MeshBase & mesh,
 }
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void connect_families(connected_elem_set_type & connected_elements,
                       const MeshBase * mesh)
 {
@@ -362,6 +363,7 @@ void connect_families(connected_elem_set_type & connected_elements,
   connect_element_families(connected_elements, connected_elements,
                            connected_elements, mesh);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/mesh/mesh_function.C b/src/mesh/mesh_function.C
index 930d8dad0a6..c1c7425ebaf 100644
--- a/src/mesh/mesh_function.C
+++ b/src/mesh/mesh_function.C
@@ -124,6 +124,7 @@ void MeshFunction::init ()
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void MeshFunction::init (const Trees::BuildType /*point_locator_build_type*/)
 {
   libmesh_deprecated();
@@ -133,6 +134,7 @@ void MeshFunction::init (const Trees::BuildType /*point_locator_build_type*/)
   // previously anyway.
   this->init();
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/reduced_basis/rb_evaluation.C b/src/reduced_basis/rb_evaluation.C
index 4ca90da8376..4e1f9690006 100644
--- a/src/reduced_basis/rb_evaluation.C
+++ b/src/reduced_basis/rb_evaluation.C
@@ -471,6 +471,7 @@ Real RBEvaluation::residual_scaling_denom(Real alpha_LB)
   return alpha_LB;
 }
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 Real RBEvaluation::eval_output_dual_norm(unsigned int n, const RBParameters & /*mu*/)
 {
   libmesh_deprecated();
@@ -478,6 +479,7 @@ Real RBEvaluation::eval_output_dual_norm(unsigned int n, const RBParameters & /*
   // Call non-deprecated version of this function, ignoring input mu
   return this->eval_output_dual_norm(n, nullptr);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 Real RBEvaluation::eval_output_dual_norm(unsigned int n,
                                          const std::vector<Number> * evaluated_thetas)
diff --git a/src/reduced_basis/rb_parametrized.C b/src/reduced_basis/rb_parametrized.C
index ac4e2d1c5f2..f756f528b28 100644
--- a/src/reduced_basis/rb_parametrized.C
+++ b/src/reduced_basis/rb_parametrized.C
@@ -126,6 +126,7 @@ unsigned int RBParametrized::get_n_discrete_params() const
     (get_discrete_parameter_values().size());
 }
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 std::set<std::string> RBParametrized::get_parameter_names() const
 {
   libmesh_deprecated();
@@ -137,6 +138,7 @@ std::set<std::string> RBParametrized::get_parameter_names() const
 
   return parameter_names;
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 bool RBParametrized::set_parameters(const RBParameters & params)
 {
diff --git a/src/systems/diff_system.C b/src/systems/diff_system.C
index ed5bf7cbe3b..d96771337cb 100644
--- a/src/systems/diff_system.C
+++ b/src/systems/diff_system.C
@@ -346,6 +346,7 @@ bool DifferentiableSystem::have_second_order_scalar_vars() const
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void DifferentiableSystem::swap_physics ( DifferentiablePhysics * & swap_physics )
 {
   // This isn't safe if users aren't very careful about memory
@@ -391,6 +392,7 @@ void DifferentiableSystem::swap_physics ( DifferentiablePhysics * & swap_physics
   // rather than just transposing
   this->disable_cache();
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/systems/equation_systems.C b/src/systems/equation_systems.C
index eb24bd95189..fe1a3558439 100644
--- a/src/systems/equation_systems.C
+++ b/src/systems/equation_systems.C
@@ -1002,12 +1002,14 @@ void EquationSystems::get_vars_active_subdomains(const std::vector<std::string>
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void EquationSystems::get_solution (std::vector<Number> & soln,
                                     std::vector<std::string> & names) const
 {
   libmesh_deprecated();
   this->build_elemental_solution_vector(soln, names);
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/src/systems/implicit_system.C b/src/systems/implicit_system.C
index 4321e279ee3..294c80869ae 100644
--- a/src/systems/implicit_system.C
+++ b/src/systems/implicit_system.C
@@ -1228,6 +1228,7 @@ std::pair<unsigned int, Real> ImplicitSystem::get_linear_solve_parameters() cons
 
 
 
+#ifdef LIBMESH_ENABLE_DEPRECATED
 void ImplicitSystem::release_linear_solver(LinearSolver<Number> *) const
 {
   // This function was originally paired with get_linear_solver()
@@ -1237,6 +1238,7 @@ void ImplicitSystem::release_linear_solver(LinearSolver<Number> *) const
   // longer needs to do any cleanup.
   libmesh_deprecated();
 }
+#endif // LIBMESH_ENABLE_DEPRECATED
 
 
 
diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h
index 488e1ecfc6c..867c654de95 100644
--- a/tests/fe/fe_test.h
+++ b/tests/fe/fe_test.h
@@ -6,6 +6,7 @@
 #include <libmesh/dof_map.h>
 #include <libmesh/elem.h>
 #include <libmesh/equation_systems.h>
+#include <libmesh/fe.h>
 #include <libmesh/fe_base.h>
 #include <libmesh/fe_interface.h>
 #include <libmesh/function_base.h>
@@ -666,9 +667,29 @@ class FETest : public FETestBase<order, family, elem_type, 1> {
 
     const FEType fe_type = this->_sys->variable_type(0);
 
+    unsigned int my_n_dofs = 0;
+
+    switch (this->_elem->dim())
+    {
+    case 0:
+      my_n_dofs = FE<0,family>::n_dofs(this->_elem, order);
+      break;
+    case 1:
+      my_n_dofs = FE<1,family>::n_dofs(this->_elem, order);
+      break;
+    case 2:
+      my_n_dofs = FE<2,family>::n_dofs(this->_elem, order);
+      break;
+    case 3:
+      my_n_dofs = FE<3,family>::n_dofs(this->_elem, order);
+      break;
+    default:
+      libmesh_error();
+    }
+
     CPPUNIT_ASSERT_EQUAL(
       FEInterface::n_shape_functions(fe_type, this->_elem),
-      this->_fe->n_shape_functions());
+      my_n_dofs);
 
     CPPUNIT_ASSERT_EQUAL(
       FEInterface::get_continuity(fe_type),
diff --git a/tests/fe/inf_fe_radial_test.C b/tests/fe/inf_fe_radial_test.C
index f8df56b0461..c6a97e5f570 100644
--- a/tests/fe/inf_fe_radial_test.C
+++ b/tests/fe/inf_fe_radial_test.C
@@ -436,9 +436,9 @@ public:
 
         const Point b = base_point(_qpoint[qp], infinite_elem) - infinite_elem->origin();
         const Point p = _qpoint[qp] - infinite_elem ->origin();
-        const Point rp = FEInterface::inverse_map(3, fe_type, infinite_elem, _qpoint[qp] );
-        const Point again_qp = FEInterface::map(3, fe_type, infinite_elem, rp);
-        //const Point rb = FEInterface::inverse_map(3, fe_type, infinite_elem, b+infinite_elem->origin());
+        const Point rp = FEMap::inverse_map(3, infinite_elem, _qpoint[qp]);
+        const Point again_qp = FEMap::map(3, infinite_elem, rp);
+        //const Point rb = FEMap::inverse_map(3, infinite_elem, b+infinite_elem->origin());
         const Real v=rp(2);
 
         // check that map() does the opposite from inverse_map
@@ -900,7 +900,7 @@ public:
 
           const std::vector<std::vector<Real> >&          phi = cfe->get_phi();
 
-          const Point node_internal = FEInterface::inverse_map(3, fe_type, some_elem, *node);
+          const Point node_internal = FEMap::inverse_map(3, some_elem, *node);
           std::vector<Point> eval_pt(1, node_internal);
           cfe->reinit(some_elem, &eval_pt);
 
@@ -932,7 +932,7 @@ public:
                 cfe = fe.get();
 
               const std::vector<std::vector<Real> >&          phi1 = cfe->get_phi();
-              eval_pt[0]=FEInterface::inverse_map(3, fe_type, elem, *node);
+              eval_pt[0]=FEMap::inverse_map(3, elem, *node);
 
               cfe->reinit(elem, &eval_pt);