Index: include/numerics/wrapped_function.h
===================================================================
--- include/numerics/wrapped_function.h	(revision 5057)
+++ include/numerics/wrapped_function.h	(working copy)
@@ -17,14 +17,15 @@
 
 
 
-#ifndef __function_base_h__
-#define __function_base_h__
+#ifndef __wrapped_function_h__
+#define __wrapped_function_h__
 
 // C++ includes
 
 
 
 // Local Includes
+#include "function_base.h"
 #include "libmesh_common.h"
 
 namespace libMesh
@@ -34,129 +35,193 @@
 
 // Forward Declarations
 template <typename T> class DenseVector;
-class FunctionBase;
 class Point;
 
 
 /**
- * This is the base class for functor-like classes.  These
- * entities are functions (in the mathematical sense) of time
- * and space, \f$ f(\mathbf{x},t) =  \mbox{\texttt{v}} \f$,
- * where \p v may be either a \p Number or a \p DenseVector<Number>.
- * Children of this base class implement different styles of
- * data retrieval for these functions.  Use the constructors
- * of the derived classes for creating new objects. The
- * required input of each derived class thwarts the effective
- * use of the commonly used \p build() member.  But afterwards
- * the virtual members allow the convenient and libMesh-common
- * usage through a \p FunctionBase*.
- *
- * @author Daniel Dreyer, 2003
  */
 
 // ------------------------------------------------------------
-// FunctionBase class definition
-class FunctionBase
+// WrappedFunction class definition
+template <typename Output=Number>
+class WrappedFunction : public FunctionBase<Output>
 {
-protected:
-
-  /**
-   * Constructor.  Optionally takes a master.
-   */
-  FunctionBase (const FunctionBase* master = NULL);
-
 public:
 
   /**
-   * Destructor.
+   * Constructor to wrap scalar-valued function pointers.
    */
-  virtual ~FunctionBase ();
+  WrappedFunction (const System &sys,
+                   Output fptr(const Point& p,
+                               const Parameters& parameters,
+                               const std::string& sys_name,
+                               const std::string& unknown_name) = NULL,
+                   const Parameters* parameters = NULL,
+                   unsigned int varnum=0)
+    : _sys(sys),
+      _parameters(parameters),
+      _fptr(fptr),
+      _varnum(varnum)
+  {
+    this->_initialized = true;
+    if (!parameters)
+      _parameters = &sys.get_equation_systems().parameters;
+  }
 
-  /**
-   * The actual initialization process.
-   */
-  virtual void init () = 0;
+  virtual void init () {}
 
-  /**
-   * Clears the function.
-   */
-  virtual void clear () = 0;
+  virtual void clear () {}
 
-
-  // ------------------------------------------------------
-  // misc
   /**
-   * @returns the scalar value at coordinate
-   * \p p and time \p time, which defaults to zero.
-   * Purely virtual, so you have to overload it.
-   * Note that this cannot be a const method, check \p MeshFunction.
+   * @returns the scalar value of variable varnum at coordinate \p p
+   * and time \p time.
    */
-  virtual Number operator() (const Point& p,
-			     const Real time = 0.) = 0;
+  virtual Output operator() (const Point& p,
+			     const Real time = 0.);
 
   /**
    * Return function for vectors.
-   * Returns in \p output the values of the data at the
-   * coordinate \p p.
-   */
-  void operator() (const Point& p,
-		   DenseVector<Number>& output);
-
-  /**
-   * Return function for vectors.
-   * Returns in \p output the values of the data at the
+   * Returns in \p output the values of all system variables at the
    * coordinate \p p and for time \p time.
-   * Purely virtual, so you have to overload it.
-   * Note that this cannot be a const method, check \p MeshFunction.
    */
   virtual void operator() (const Point& p,
 			   const Real time,
-			   DenseVector<Number>& output) = 0;
+			   DenseVector<Output>& output);
 
   /**
-   * @returns \p true when this object is properly initialized
-   * and ready for use, \p false otherwise.
+   * @returns the vector component \p i at coordinate
+   * \p p and time \p time.
    */
-  bool initialized () const;
+  virtual Output operator() (const unsigned int i,
+                             const Point& p,
+			     const Real time=0.);
 
-
 protected:
 
-  /**
-   * Const pointer to our master, initialized to \p NULL.
-   * There may be cases where multiple functions are required,
-   * but to save memory, one master handles some centralized
-   * data.
-   */
-  const FunctionBase* _master;
+  const System& _sys;
 
-  /**
-   * When \p init() was called so that everything is ready
-   * for calls to \p operator() (...), then this \p bool is true.
-   */
-  bool _initialized;
+  const Parameters* _parameters;
 
+  Output (*_fptr)(const Point& p,
+                  const Parameters& parameters,
+                  const std::string& sys_name,
+                  const std::string& unknown_name);
+
+  unsigned int _varnum;
 };
 
 
 // ------------------------------------------------------------
-// FunctionBase inline methods
+// WrappedFunction inline methods
+
+
+template <typename Output>
 inline
-bool FunctionBase::initialized() const
+Output WrappedFunction<Output>::operator() (const Point& p,
+			                    const Real time)
 {
-  return (this->_initialized);
+  libmesh_assert(_fptr);
+  libmesh_assert(_parameters);
+  return _fptr(p,
+               *_parameters,
+               _sys.name(),
+               _sys.variable_name(_varnum));
 }
 
+/**
+ * Return function for vectors.
+ * Returns in \p output the values of all system variables at the
+ * coordinate \p p and for time \p time.
+ */
+template <typename Output>
+inline
+void WrappedFunction<Output>::operator() (const Point& p,
+			                  const Real time,
+			                  DenseVector<Output>& output)
+{
+  libmesh_assert(_fptr);
+  libmesh_assert(_parameters);
 
+  // We fill each entry of output with a single scalar component of
+  // the data in our System
+  const unsigned int size = output.size();
+  libmesh_assert(size == _sys.n_components());
 
+  // Loop over variables, then over each component in
+  // vector-valued variables, evaluating each.
+  const unsigned int n_vars = _sys.n_vars();
+  for (unsigned int v = 0; v != n_vars; ++v)
+    {
+      const unsigned int n_components =
+        _sys.variable(v).n_components();
+      if (n_components == 1)
+        output(_sys.variable_scalar_number(v,0)) =
+	   _fptr(p, *_parameters, _sys.name(), _sys.variable_name(v));
+      else
+        {
+          // Right now our only non-scalar variable type is the
+          // SCALAR variables.  The irony is priceless.
+          libmesh_assert(_sys.variable(v).type().family == SCALAR);
+
+	  // We pass the point (j,0,0) to an old-style fptr function
+	  // pointer to distinguish the different scalars within the
+	  // SCALAR variable.
+          for (unsigned int j=1; j != n_components; ++j)
+            output(_sys.variable_scalar_number(v,j)) =
+	       _fptr(Point(j,0,0), *_parameters,
+                     _sys.name(), _sys.variable_name(v));
+        }
+    }
+}
+
+
+/**
+ * @returns the vector component \p i at coordinate
+ * \p p and time \p time.
+ */
+template <typename Output>
 inline
-void FunctionBase::operator() (const Point& p,
-			       DenseVector<Number>& output)
+Output WrappedFunction<Output>::operator() (const unsigned int i,
+                                            const Point& p,
+			                    const Real time)
 {
-  // Call the time-dependent function with t=0.
-  this->operator()(p, 0., output);
+  libmesh_assert(_fptr);
+  libmesh_assert(_parameters);
+
+  // Loop over variables, then over each component in
+  // vector-valued variables.
+  const unsigned int n_vars = _sys.n_vars();
+  for (unsigned int v = 0; v != n_vars; ++v)
+    {
+      const unsigned int n_components =
+        _sys.variable(v).n_components();
+      if (n_components == 1 &&
+          i == _sys.variable_scalar_number(v,0))
+        return _fptr(p, *_parameters, _sys.name(), _sys.variable_name(v));
+      else if (i >= _sys.variable_scalar_number(v,0) &&
+               i <= _sys.variable_scalar_number(v,n_components-1))
+        {
+          // Right now our only non-scalar variable type is the
+          // SCALAR variables.  The irony is priceless.
+          libmesh_assert(_sys.variable(i).type().family == SCALAR);
+
+	  // We pass the point (j,0,0) to an old-style fptr function
+	  // pointer to distinguish the different scalars within the
+	  // SCALAR variable.
+          for (unsigned int j=1; j != n_components; ++j)
+            if (i == _sys.variable_scalar_number(v,j))
+              return _fptr(Point(j,0,0), *_parameters,
+			   _sys.name(), _sys.variable_name(j));
+        }
+    }
+  libMesh::err << "Component index " << i <<
+		  " not found in system " << _sys.name() << std::endl;
+  libmesh_error();
+  return Output();
 }
 
+
+
 } // namespace libMesh
 
-#endif
+#endif // __wrapped_function__
Index: include/numerics/analytic_function.h
===================================================================
--- include/numerics/analytic_function.h	(revision 5063)
+++ include/numerics/analytic_function.h	(working copy)
@@ -46,7 +46,8 @@
 
 // ------------------------------------------------------------
 // AnalyticFunction class definition
-class AnalyticFunction : public FunctionBase
+template <typename Output=Number>
+class AnalyticFunction : public FunctionBase<Output>
 {
 public:
 
@@ -54,14 +55,14 @@
    * Constructor.  Takes a function pointer for scalar
    * return values.
    */
-  AnalyticFunction (Number fptr(const Point& p,
+  AnalyticFunction (Output fptr(const Point& p,
 				const Real time));
 
   /**
    * Constructor.  Takes a function pointer for
    * vector valued functions.
    */
-  AnalyticFunction (void fptr(DenseVector<Number>& output,
+  AnalyticFunction (void fptr(DenseVector<Output>& output,
 			      const Point& p,
 			      const Real time));
   /**
@@ -75,13 +76,13 @@
    * the boundary values when an analytical expression
    * is available.
    */
-  Number (* _number_fptr) (const Point& p,
+  Output (* _number_fptr) (const Point& p,
 			   const Real time);
 
   /**
    * Pointer to user-provided vector valued function.
    */
-  void (* _vector_fptr) (DenseVector<Number>& output,
+  void (* _vector_fptr) (DenseVector<Output>& output,
 			 const Point& p,
 			 const Real time);
 
@@ -99,7 +100,7 @@
    * @returns the value at point \p p and time
    * \p time, which defaults to zero.
    */
-  Number operator() (const Point& p,
+  Output operator() (const Point& p,
 		     const Real time=0.);
 
   /**
@@ -108,7 +109,7 @@
    */
   void operator() (const Point& p,
 		   const Real time,
-		   DenseVector<Number>& output);
+		   DenseVector<Output>& output);
 
 };
 
@@ -116,9 +117,10 @@
 
 // ------------------------------------------------------------
 // AnalyticFunction inline methods
+template <typename Output>
 inline
-Number AnalyticFunction::operator() (const Point& p,
-				     const Real time)
+Output AnalyticFunction<Output>::operator() (const Point& p,
+				             const Real time)
 {
   libmesh_assert (this->initialized());
   return (this->_number_fptr(p, time));
@@ -126,10 +128,11 @@
 
 
 
+template <typename Output>
 inline
-void AnalyticFunction::operator() (const Point& p,
-				   const Real time,
-				   DenseVector<Number>& output)
+void AnalyticFunction<Output>::operator() (const Point& p,
+				           const Real time,
+				           DenseVector<Output>& output)
 {
   libmesh_assert (this->initialized());
   this->_vector_fptr(output, p, time);
@@ -137,6 +140,66 @@
 }
 
 
+
+template <typename Output>
+AnalyticFunction<Output>::AnalyticFunction (Output fptr(const Point& p,
+						        const Real time)) :
+  FunctionBase<Output> (),
+  _number_fptr (fptr)
+{
+  libmesh_assert (fptr != NULL);
+  this->_initialized = true;
+}
+
+
+
+template <typename Output>
+inline
+AnalyticFunction<Output>::AnalyticFunction (void fptr(DenseVector<Output>& output,
+					              const Point& p,
+					              const Real time)) :
+  FunctionBase<Output> (),
+  _vector_fptr (fptr)
+{
+  libmesh_assert (fptr != NULL);
+  this->_initialized = true;
+}
+
+
+
+template <typename Output>
+inline
+AnalyticFunction<Output>::~AnalyticFunction ()
+{
+}
+
+
+
+template <typename Output>
+void AnalyticFunction<Output>::init ()
+{
+  // dumb double-test
+  libmesh_assert ((_number_fptr != NULL) || (_vector_fptr != NULL));
+
+  // definitely ready
+  this->_initialized = true;
+}
+
+
+
+template <typename Output>
+inline
+void AnalyticFunction<Output>::clear ()
+{
+  // We probably need a method to reset these later...
+  _number_fptr = NULL;
+  _vector_fptr = NULL;
+
+  // definitely not ready
+  this->_initialized = false;
+}
+
+
 } // namespace libMesh
 
 
Index: include/numerics/function_base.h
===================================================================
--- include/numerics/function_base.h	(revision 5063)
+++ include/numerics/function_base.h	(working copy)
@@ -34,7 +34,6 @@
 
 // Forward Declarations
 template <typename T> class DenseVector;
-class FunctionBase;
 class Point;
 
 
@@ -56,6 +55,7 @@
 
 // ------------------------------------------------------------
 // FunctionBase class definition
+template <typename Output=Number>
 class FunctionBase
 {
 protected:
@@ -91,7 +91,7 @@
    * Purely virtual, so you have to overload it.
    * Note that this cannot be a const method, check \p MeshFunction.
    */
-  virtual Number operator() (const Point& p,
+  virtual Output operator() (const Point& p,
 			     const Real time = 0.) = 0;
 
   /**
@@ -100,7 +100,7 @@
    * coordinate \p p.
    */
   void operator() (const Point& p,
-		   DenseVector<Number>& output);
+		   DenseVector<Output>& output);
 
   /**
    * Return function for vectors.
@@ -111,9 +111,24 @@
    */
   virtual void operator() (const Point& p,
 			   const Real time,
-			   DenseVector<Number>& output) = 0;
+			   DenseVector<Output>& output) = 0;
 
   /**
+   * @returns the vector component \p i at coordinate
+   * \p p and time \p time.
+   * Subclasses aren't required to overload this, since the default
+   * implementation is based on the full vector evaluation, which is
+   * often correct.
+   * Subclasses are recommended to overload this, since the default
+   * implementation is based on a vector evaluation, which is usually
+   * unnecessarily inefficient.
+   */
+  virtual Output operator() (const unsigned int i,
+                             const Point& p,
+			     const Real time=0.);
+
+
+  /**
    * @returns \p true when this object is properly initialized
    * and ready for use, \p false otherwise.
    */
@@ -141,18 +156,50 @@
 
 // ------------------------------------------------------------
 // FunctionBase inline methods
+
+template<typename Output>
 inline
-bool FunctionBase::initialized() const
+FunctionBase<Output>::FunctionBase (const FunctionBase* master) :
+  _master             (master),
+  _initialized        (false)
 {
+}
+
+
+
+template<typename Output>
+inline
+FunctionBase<Output>::~FunctionBase ()
+{
+}
+
+
+
+template <typename Output>
+inline
+bool FunctionBase<Output>::initialized() const
+{
   return (this->_initialized);
 }
 
 
 
+template <typename Output>
 inline
-void FunctionBase::operator() (const Point& p,
-			       DenseVector<Number>& output)
+Output FunctionBase<Output>::operator() (const unsigned int i,
+                                         const Point& p,
+			                 const Real time)
 {
+  DenseVector<Output> out(i+1);
+  (*this)(p, time, out);
+  return out(i);
+}
+
+template <typename Output>
+inline
+void FunctionBase<Output>::operator() (const Point& p,
+			               DenseVector<Output>& output)
+{
   // Call the time-dependent function with t=0.
   this->operator()(p, 0., output);
 }
Index: include/numerics/parsed_function.h
===================================================================
--- include/numerics/parsed_function.h	(revision 5063)
+++ include/numerics/parsed_function.h	(working copy)
@@ -13,7 +13,8 @@
 
 #include "fparser.hh"
 
-class ParsedFunction : public FunctionBase
+template <typename Output=Number>
+class ParsedFunction : public FunctionBase<Output>
 {
 public:
   ParsedFunction (const std::string& expression)
@@ -42,7 +43,7 @@
         start = (end == std::string::npos) ?
                 std::string::npos : end + 1;
 
-        FunctionParser fp;
+        FunctionParserBase<Output> fp;
         fp.AddConstant("pi", std::acos(Real(-1)));
         fp.AddConstant("e", std::exp(Real(1)));
         fp.Parse(subexpression, variables);
@@ -50,7 +51,7 @@
       }
     }
 
-  virtual Number operator() (const Point& p,
+  virtual Output operator() (const Point& p,
                              const Real time = 0)
     {
       _spacetime[0] = p(0);
@@ -66,7 +67,7 @@
 
   virtual void operator() (const Point& p,
                            const Real time,
-                           DenseVector<Number>& output)
+                           DenseVector<Output>& output)
     {
       _spacetime[0] = p(0);
 #if LIBMESH_DIM > 1
@@ -88,13 +89,14 @@
   virtual void clear() {}
 
 private:
-  Number _spacetime[LIBMESH_DIM+1];
-  std::vector<FunctionParserBase<Number> > parsers;
+  Output _spacetime[LIBMESH_DIM+1];
+  std::vector<FunctionParserBase<Output> > parsers;
 };
 
 #else
 
-class ParsedFunction : public FunctionBase
+template <typename Output>
+class ParsedFunction : public FunctionBase<Output>
 {
 public:
   ParsedFunction (std::string expression)
@@ -102,13 +104,16 @@
       libmesh_not_implemented();
     }
 
-  virtual Number operator() (const Point&,
+  virtual Output operator() (const Point&,
                              const Real time = 0)
     { return 0.; }
 
   virtual void operator() (const Point&,
                            const Real time,
-                           DenseVector<Number>& output) {}
+                           DenseVector<Output>& output) {}
+
+  virtual void init() {}
+  virtual void clear() {}
 };
 
 #endif // LIBMESH_HAVE_FPARSER
Index: include/numerics/zero_function.h
===================================================================
--- include/numerics/zero_function.h	(revision 5063)
+++ include/numerics/zero_function.h	(working copy)
@@ -8,9 +8,10 @@
 #include "function_base.h"
 #include "point.h"
 
-class ZeroFunction : public FunctionBase
+template <typename Output=Number>
+class ZeroFunction : public FunctionBase<Output>
 {
-  virtual Number operator() (const Point&,
+  virtual Output operator() (const Point&,
                              const Real = 0)
     {
       return 0;
@@ -18,7 +19,7 @@
 
   virtual void operator() (const Point&,
                            const Real,
-                           DenseVector<Number>& output)
+                           DenseVector<Output>& output)
     {
       unsigned int size = output.size();
       for (unsigned int i=0; i != size; ++i)
Index: include/mesh/mesh_function.h
===================================================================
--- include/mesh/mesh_function.h	(revision 5063)
+++ include/mesh/mesh_function.h	(working copy)
@@ -53,7 +53,7 @@
 
 // ------------------------------------------------------------
 // MeshFunction class definition
-class MeshFunction : public FunctionBase
+class MeshFunction : public FunctionBase<Number>
 {
 public:
 
@@ -69,7 +69,7 @@
 		const NumericVector<Number>& vec,
 		const DofMap& dof_map,
 		const std::vector<unsigned int>& vars,
-		const FunctionBase* master=NULL);
+		const FunctionBase<>* master=NULL);
 
   /**
    * Constructor for mesh based functions with a number
@@ -83,7 +83,7 @@
 		const NumericVector<Number>& vec,
 		const DofMap& dof_map,
 		const unsigned int var,
-		const FunctionBase* master=NULL);
+		const FunctionBase<>* master=NULL);
 
   /**
    * Destructor.
Index: include/systems/system.h
===================================================================
--- include/systems/system.h	(revision 5063)
+++ include/systems/system.h	(working copy)
@@ -48,6 +48,7 @@
 class MeshBase;
 class Xdr;
 class DofMap;
+template <typename Output> class FunctionBase;
 class Parameters;
 class ParameterVector;
 class Point;
@@ -443,7 +444,11 @@
   virtual std::string system_type () const { return "Basic"; }
 
   /**
-   * Projects the continuous functions onto the current solution.
+   * Projects arbitrary functions onto the current solution.
+   * The function value \p fptr and its gradient \p gptr are
+   * represented by function pointers.
+   * A gradient \p gptr is only required/used for projecting onto
+   * finite element spaces with continuous derivatives.
    */
   void project_solution (Number fptr(const Point& p,
 				     const Parameters& parameters,
@@ -456,8 +461,26 @@
 			 Parameters& parameters) const;
 
   /**
-   * Projects the continuous functions onto the current mesh.
+   * Projects arbitrary functions onto the current solution.
+   * The function value \p f and its gradient \p g are
+   * represented by functors.
+   * A gradient \p g is only required/used for projecting onto finite
+   * element spaces with continuous derivatives.
+   * If non-default \p Parameters are to be used, they can be provided
+   * in the \p parameters argument.
    */
+  void project_solution (FunctionBase<Number> *f,
+                         FunctionBase<Gradient> *g = NULL,
+                         Parameters* parameters = NULL) const;
+
+  /**
+   * Projects arbitrary functions onto a vector of degree of freedom
+   * values for the current system.
+   * The function value \p fptr and its gradient \p gptr are
+   * represented by function pointers.
+   * A gradient \p gptr is only required/used for projecting onto
+   * finite element spaces with continuous derivatives.
+   */
   void project_vector (Number fptr(const Point& p,
 				   const Parameters& parameters,
                                    const std::string& sys_name,
@@ -470,6 +493,21 @@
 		       NumericVector<Number>& new_vector) const;
 
   /**
+   * Projects arbitrary functions onto a vector of degree of freedom
+   * values for the current system.
+   * The function value \p f and its gradient \p g are
+   * represented by functors.
+   * A gradient \p g is only required/used for projecting onto finite
+   * element spaces with continuous derivatives.
+   * If non-default \p Parameters are to be used, they can be provided
+   * in the \p parameters argument.
+   */
+  void project_vector (NumericVector<Number>& new_vector,
+                       FunctionBase<Number> *f,
+                       FunctionBase<Gradient> *g = NULL,
+		       Parameters* parameters = NULL) const;
+
+  /**
    * @returns the system number.
    */
   unsigned int number () const;
@@ -796,6 +834,13 @@
   unsigned int n_vars() const;
 
   /**
+   * @returns the total number of scalar components in the system's
+   * variables.  This will equal \p n_vars() in the case of all
+   * scalar-valued variables.
+   */
+  unsigned int n_components() const;
+
+  /**
    * @returns the number of degrees of freedom in the system
    */
   unsigned int n_dofs() const;
@@ -858,6 +903,33 @@
   unsigned short int variable_number (const std::string& var) const;
 
   /**
+   * @returns an index, starting from 0 for the first component of the
+   * first variable, and incrementing for each component of each
+   * (potentially vector-valued) variable in the system in order.
+   * For systems with only scalar-valued variables, this will be the
+   * same as \p variable_number(var)
+   *
+   * Irony: currently our only non-scalar-valued variable type is
+   * SCALAR.
+   */
+  unsigned short int variable_scalar_number (const std::string& var,
+                                             unsigned short component) const;
+
+  /**
+   * @returns an index, starting from 0 for the first component of the
+   * first variable, and incrementing for each component of each
+   * (potentially vector-valued) variable in the system in order.
+   * For systems with only scalar-valued variables, this will be the
+   * same as \p var_num
+   *
+   * Irony: currently our only non-scalar-valued variable type is
+   * SCALAR.
+   */
+  unsigned short int variable_scalar_number (unsigned short var_num,
+                                             unsigned short component) const;
+
+
+  /**
    * @returns the finite element type variable number \p i.
    */
   const FEType & variable_type (const unsigned int i) const;
@@ -1453,6 +1525,12 @@
   std::map<std::string, unsigned short int> _variable_numbers;
 
   /**
+   * The variable scalar numbers for the first component of each
+   * variable.
+   */
+  std::vector<unsigned short> _variable_first_scalar;
+
+  /**
    * Flag stating if the system is active or not.
    */
   bool _active;
@@ -1581,34 +1659,20 @@
   private:
     const System                &system;
 
-    Number (* fptr)(const Point& p,
-		    const Parameters& parameters,
-		    const std::string& sys_name,
-		    const std::string& unknown_name);
-
-    Gradient (* gptr)(const Point& p,
-		      const Parameters& parameters,
-		      const std::string& sys_name,
-		      const std::string& unknown_name);
-
-    Parameters &parameters;
+    FunctionBase<Number> *f;
+    FunctionBase<Gradient> *g;
+    const Parameters &parameters;
     NumericVector<Number>       &new_vector;
 
   public:
     ProjectSolution (const System &system_in,
-		     Number fptr_in(const Point& p,
-				    const Parameters& parameters,
-				    const std::string& sys_name,
-				    const std::string& unknown_name),
-		     Gradient gptr_in(const Point& p,
-				      const Parameters& parameters,
-				      const std::string& sys_name,
-				      const std::string& unknown_name),
-		     Parameters &parameters_in,
+		     FunctionBase<Number>* f_in,
+		     FunctionBase<Gradient>* g_in,
+		     const Parameters &parameters_in,
 		     NumericVector<Number> &new_v_in) :
     system(system_in),
-    fptr(fptr_in),
-    gptr(gptr_in),
+    f(f_in),
+    g(g_in),
     parameters(parameters_in),
     new_vector(new_v_in)
     {}
@@ -1710,6 +1774,14 @@
 
 
 inline
+unsigned int System::n_components() const
+{
+  return _variable_first_scalar[n_vars()];
+}
+
+
+
+inline
 const Variable & System::variable (const unsigned int i) const
 {
   libmesh_assert (i < _variables.size());
@@ -1730,6 +1802,26 @@
 
 
 inline
+unsigned short int 
+System::variable_scalar_number (const std::string& var,
+                                unsigned short component) const
+{
+  return variable_scalar_number(this->variable_number(var), component);
+}
+
+
+
+inline
+unsigned short int
+System::variable_scalar_number (unsigned short var_num, 
+                                unsigned short component) const
+{
+  return _variable_first_scalar[var_num] + component;
+}
+
+
+
+inline
 const FEType & System::variable_type (const unsigned int i) const
 {
   libmesh_assert (i < _variables.size());
Index: include/base/variable.h
===================================================================
--- include/base/variable.h	(revision 5063)
+++ include/base/variable.h	(working copy)
@@ -87,6 +87,12 @@
   { return _type; }
 
   /**
+   * The number of components of this variable.
+   */
+  unsigned int n_components() const
+  { return type().family == SCALAR ? _type.order : 1; }
+
+  /**
    * \p returns \p true if this variable is active on subdomain \p sid,
    * \p false otherwise.  Note that we interperet the special case of an
    * empty \p _active_subdomains container as active everywhere, i.e.
Index: src/numerics/mesh_function.C
===================================================================
--- src/numerics/mesh_function.C	(revision 5063)
+++ src/numerics/mesh_function.C	(working copy)
@@ -43,8 +43,8 @@
 			    const NumericVector<Number>& vec,
 			    const DofMap& dof_map,
 			    const std::vector<unsigned int>& vars,
-			    const FunctionBase* master) :
-  FunctionBase   (master),
+			    const FunctionBase<>* master) :
+  FunctionBase<> (master),
   _eqn_systems   (eqn_systems),
   _vector        (vec),
   _dof_map       (dof_map),
@@ -61,8 +61,8 @@
 			    const NumericVector<Number>& vec,
 			    const DofMap& dof_map,
 			    const unsigned int var,
-			    const FunctionBase* master) :
-  FunctionBase   (master),
+			    const FunctionBase<>* master) :
+  FunctionBase<> (master),
   _eqn_systems   (eqn_systems),
   _vector        (vec),
   _dof_map       (dof_map),
Index: src/numerics/analytic_function.C
===================================================================
--- src/numerics/analytic_function.C	(revision 5063)
+++ src/numerics/analytic_function.C	(working copy)
@@ -25,64 +25,4 @@
 
 namespace libMesh
 {
-
-
-
-//------------------------------------------------------------------
-// AnalyticFunction methods
-AnalyticFunction::AnalyticFunction (Number fptr(const Point& p,
-						const Real time)) :
-  FunctionBase (),
-  _number_fptr (fptr)
-{
-  libmesh_assert (fptr != NULL);
-  this->_initialized = true;
-}
-
-
-
-
-AnalyticFunction::AnalyticFunction (void fptr(DenseVector<Number>& output,
-					      const Point& p,
-					      const Real time)) :
-  FunctionBase (),
-  _vector_fptr (fptr)
-{
-  libmesh_assert (fptr != NULL);
-  this->_initialized = true;
-}
-
-
-
-
-AnalyticFunction::~AnalyticFunction ()
-{
-}
-
-
-
-
-
-void AnalyticFunction::init ()
-{
-  // dumb double-test
-  libmesh_assert ((_number_fptr != NULL) || (_vector_fptr != NULL));
-
-  // definitely ready
-  this->_initialized = true;
-}
-
-
-
-void AnalyticFunction::clear ()
-{
-  // We probably need a method to reset these later...
-  _number_fptr = NULL;
-  _vector_fptr = NULL;
-
-  // definitely not ready
-  this->_initialized = false;
-}
-
-
 } // namespace libMesh
Index: src/numerics/function_base.C
===================================================================
--- src/numerics/function_base.C	(revision 5063)
+++ src/numerics/function_base.C	(working copy)
@@ -28,21 +28,6 @@
 
 
 
-
-//------------------------------------------------------------------
-// FunctionBase methods
-FunctionBase::FunctionBase (const FunctionBase* master) :
-  _master             (master),
-  _initialized        (false)
-{
-}
-
-
-
-FunctionBase::~FunctionBase ()
-{
-}
-
 } // namespace libMesh
 
 
Index: src/systems/system_projection.C
===================================================================
--- src/systems/system_projection.C	(revision 5063)
+++ src/systems/system_projection.C	(working copy)
@@ -28,12 +28,14 @@
 #include "fe_interface.h"
 #include "numeric_vector.h"
 #include "libmesh_logging.h"
+#include "equation_systems.h"
 
 #include "dense_matrix.h"
 #include "fe_base.h"
 #include "quadrature_gauss.h"
 #include "dense_vector.h"
 #include "threads.h"
+#include "wrapped_function.h"
 
 namespace libMesh
 {
@@ -232,7 +234,7 @@
 
 
 /**
- * This method projects an analytic function onto the solution via L2
+ * This method projects an arbitrary function onto the solution via L2
  * projections and nodal interpolations on each element.
  */
 void System::project_solution (Number fptr(const Point& p,
@@ -245,13 +247,29 @@
                                              const std::string& unknown_name),
                                Parameters& parameters) const
 {
-  this->project_vector(fptr, gptr, parameters, *solution);
+  WrappedFunction<Number> f(*this, fptr, &parameters);
+  WrappedFunction<Gradient> g(*this, gptr, &parameters);
+  this->project_solution(&f, &g);
+}
 
+
+/**
+ * This method projects an analytic function onto the solution via L2
+ * projections and nodal interpolations on each element.
+ */
+void System::project_solution (FunctionBase<Number> *f,
+                               FunctionBase<Gradient> *g,
+                               Parameters* parameters) const
+{
+  this->project_vector(*solution, f, g, parameters);
+
   solution->localize(*current_local_solution);
 }
 
 
 
+
+
 /**
  * This method projects an analytic function via L2 projections and
  * nodal interpolations on each element.
@@ -267,23 +285,43 @@
                              Parameters& parameters,
                              NumericVector<Number>& new_vector) const
 {
+  WrappedFunction<Number> f(*this, fptr, &parameters);
+  WrappedFunction<Gradient> g(*this, gptr, &parameters);
+  this->project_vector(new_vector, &f, &g);
+}
+
+/**
+ * This method projects an analytic function via L2 projections and
+ * nodal interpolations on each element.
+ */
+void System::project_vector (NumericVector<Number>& new_vector,
+                             FunctionBase<Number> *f,
+                             FunctionBase<Gradient> *g,
+                             Parameters *parameters) const
+{
   START_LOG ("project_vector()", "System");
 
-  Threads::parallel_for (ConstElemRange (this->get_mesh().active_local_elements_begin(),
-					 this->get_mesh().active_local_elements_end(),
-					 1000),
-			 ProjectSolution(*this,
-					 fptr,
-					 gptr,
-					 parameters,
-					 new_vector)
-			 );
+  Threads::parallel_for
+    (ConstElemRange (this->get_mesh().active_local_elements_begin(),
+		     this->get_mesh().active_local_elements_end(),
+		     1000),
+     ProjectSolution(*this, f, g,
+		     parameters ?
+		       *parameters :
+                       this->get_equation_systems().parameters,
+		     new_vector)
+    );
 
   // Also, load values into the SCALAR dofs
   // Note: We assume that all SCALAR dofs are on the
   // processor with highest ID
   if(libMesh::processor_id() == (libMesh::n_processors()-1))
   {
+    // We get different scalars as different
+    // components from a new-style f functor.
+    DenseVector<Number> fout(this->n_components());
+    (*f) (Point(), 0, fout);
+
     const DofMap& dof_map = this->get_dof_map();
     for (unsigned int var=0; var<this->n_vars(); var++)
       if(this->variable(var).type().family == SCALAR)
@@ -294,16 +332,10 @@
 
           for (unsigned int i=0; i<n_SCALAR_dofs; i++)
           {
-            const unsigned int index = SCALAR_indices[i];
-
-            // We pass the point (i,0,0) to the fptr to distinguish
-            // the different scalars within the SCALAR variable
-            Point p_i(i,0,0);
-            new_vector.set( index, fptr(p_i,
-                                        parameters,
-                                        this->name(),
-                                        this->variable_name(var))
-                          );
+            const unsigned int global_index = SCALAR_indices[i];
+            const unsigned int component_index =
+              this->variable_scalar_number(var,i);
+            new_vector.set(global_index, fout(component_index));
           }
         }
   }
@@ -322,6 +354,7 @@
 void System::ProjectVector::operator()(const ConstElemRange &) const
 {
   libmesh_error();
+}
 #else
 void System::ProjectVector::operator()(const ConstElemRange &range) const
 {
@@ -735,9 +768,8 @@
     } // end variables loop
 
   STOP_LOG ("operator()","ProjectVector");
-
-#endif // LIBMESH_ENABLE_AMR
 }
+#endif // LIBMESH_ENABLE_AMR
 
 
 
@@ -765,6 +797,7 @@
 void System::BuildProjectionList::operator()(const ConstElemRange &)
 {
   libmesh_error();
+}
 #else
 void System::BuildProjectionList::operator()(const ConstElemRange &range)
 {
@@ -833,8 +866,8 @@
         if (di[i] < first_old_dof || di[i] >= end_old_dof)
           this->send_list.push_back(di[i]);
     }  // end elem loop
-#endif // LIBMESH_ENABLE_AMR
 }
+#endif // LIBMESH_ENABLE_AMR
 
 
 
@@ -850,11 +883,11 @@
 void System::ProjectSolution::operator()(const ConstElemRange &range) const
 {
   // We need data to project
-  libmesh_assert(fptr);
+  libmesh_assert(f);
 
   /**
    * This method projects an analytic solution to the current
-   * mesh.  The input function \p fptr gives the analytic solution,
+   * mesh.  The input function \p f gives the analytic solution,
    * while the \p new_vector (which should already be correctly sized)
    * gives the solution (to be computed) on the current mesh.
    */
@@ -888,6 +921,9 @@
       if (fe_type.family == SCALAR)
         continue;
 
+      const unsigned int var_component =
+        system.variable_scalar_number(var, 0);
+
       // Get FE objects of the appropriate type
       AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
 
@@ -909,7 +945,7 @@
       if (cont == C_ONE)
         {
           // We'll need gradient data for a C1 projection
-          libmesh_assert(gptr);
+          libmesh_assert(g);
 
           const std::vector<std::vector<RealGradient> >&
             ref_dphi = fe->get_dphi();
@@ -989,28 +1025,19 @@
               else if (cont == C_ZERO)
                 {
                   libmesh_assert(nc == 1);
-		  Ue(current_dof) = fptr(elem->point(n),
-                                         parameters,
-                                         system.name(),
-                                         system.variable_name(var));
+		  Ue(current_dof) = (*f)(var_component,elem->point(n));
                   dof_is_fixed[current_dof] = true;
                   current_dof++;
                 }
               // The hermite element vertex shape functions are weird
               else if (fe_type.family == HERMITE)
                 {
-                  Ue(current_dof) = fptr(elem->point(n),
-                                         parameters,
-                                         system.name(),
-                                         system.variable_name(var));
+                  Ue(current_dof) = (*f)(var_component,elem->point(n));
                   dof_is_fixed[current_dof] = true;
                   current_dof++;
-                  Gradient g = gptr(elem->point(n),
-                                    parameters,
-                                    system.name(),
-                                    system.variable_name(var));
+                  Gradient grad = (*g)(var_component,elem->point(n));
                   // x derivative
-                  Ue(current_dof) = g(0);
+                  Ue(current_dof) = grad(0);
                   dof_is_fixed[current_dof] = true;
                   current_dof++;
                   if (dim > 1)
@@ -1020,16 +1047,10 @@
                             nxplus = elem->point(n);
                       nxminus(0) -= TOLERANCE;
                       nxplus(0) += TOLERANCE;
-                      Gradient gxminus = gptr(nxminus,
-                                              parameters,
-                                              system.name(),
-                                              system.variable_name(var));
-                      Gradient gxplus = gptr(nxminus,
-                                             parameters,
-                                             system.name(),
-                                             system.variable_name(var));
+                      Gradient gxminus = (*g)(var_component,nxminus);
+                      Gradient gxplus = (*g)(var_component,nxplus);
                       // y derivative
-                      Ue(current_dof) = g(1);
+                      Ue(current_dof) = grad(1);
                       dof_is_fixed[current_dof] = true;
                       current_dof++;
                       // xy derivative
@@ -1041,7 +1062,7 @@
                       if (dim > 2)
                         {
                           // z derivative
-                          Ue(current_dof) = g(2);
+                          Ue(current_dof) = grad(2);
                           dof_is_fixed[current_dof] = true;
                           current_dof++;
                           // xz derivative
@@ -1054,14 +1075,8 @@
                                 nyplus = elem->point(n);
                           nyminus(1) -= TOLERANCE;
                           nyplus(1) += TOLERANCE;
-                          Gradient gyminus = gptr(nyminus,
-                                                  parameters,
-                                                  system.name(),
-                                                  system.variable_name(var));
-                          Gradient gyplus = gptr(nyminus,
-                                                 parameters,
-                                                 system.name(),
-                                                 system.variable_name(var));
+                          Gradient gyminus = (*g)(var_component,nyminus);
+                          Gradient gyplus = (*g)(var_component,nyplus);
                           // xz derivative
                           Ue(current_dof) = (gyplus(2) - gyminus(2))
                                             / 2. / TOLERANCE;
@@ -1080,22 +1095,10 @@
                           nxpym(1) -= TOLERANCE;
                           nxpyp(0) += TOLERANCE;
                           nxpyp(1) += TOLERANCE;
-                          Gradient gxmym = gptr(nxmym,
-                                                parameters,
-                                                system.name(),
-                                                system.variable_name(var));
-                          Gradient gxmyp = gptr(nxmyp,
-                                                parameters,
-                                                system.name(),
-                                                system.variable_name(var));
-                          Gradient gxpym = gptr(nxpym,
-                                                parameters,
-                                                system.name(),
-                                                system.variable_name(var));
-                          Gradient gxpyp = gptr(nxpyp,
-                                                parameters,
-                                                system.name(),
-                                                system.variable_name(var));
+                          Gradient gxmym = (*g)(var_component,nxmym);
+                          Gradient gxmyp = (*g)(var_component,nxmyp);
+                          Gradient gxpym = (*g)(var_component,nxpym);
+                          Gradient gxpyp = (*g)(var_component,nxpyp);
                           Number gxzplus = (gxpyp(2) - gxmyp(2))
                                          / 2. / TOLERANCE;
                           Number gxzminus = (gxpym(2) - gxmym(2))
@@ -1114,19 +1117,13 @@
               else if (cont == C_ONE)
                 {
                   libmesh_assert(nc == 1 + dim);
-		  Ue(current_dof) = fptr(elem->point(n),
-                                         parameters,
-                                         system.name(),
-                                         system.variable_name(var));
+		  Ue(current_dof) = (*f)(var_component,elem->point(n));
                   dof_is_fixed[current_dof] = true;
                   current_dof++;
-                  Gradient g = gptr(elem->point(n),
-                                    parameters,
-                                    system.name(),
-                                    system.variable_name(var));
+                  Gradient grad = (*g)(var_component,elem->point(n));
                   for (unsigned int i=0; i!= dim; ++i)
                     {
-		      Ue(current_dof) = g(i);
+		      Ue(current_dof) = grad(i);
                       dof_is_fixed[current_dof] = true;
                       current_dof++;
                     }
@@ -1167,16 +1164,11 @@
 	        for (unsigned int qp=0; qp<n_qp; qp++)
 	          {
 	            // solution at the quadrature point
-	            Number fineval = fptr(xyz_values[qp],
-                                          parameters,
-                                          system.name(),
-                                          system.variable_name(var));
+	            Number fineval = (*f)(var_component,xyz_values[qp]);
 	            // solution grad at the quadrature point
 	            Gradient finegrad;
                     if (cont == C_ONE)
-                      finegrad = gptr(xyz_values[qp], parameters,
-                                      system.name(),
-				      system.variable_name(var));
+                      finegrad = (*g)(var_component,xyz_values[qp]);
 
                     // Form edge projection matrix
                     for (unsigned int sidei=0, freei=0;
@@ -1263,16 +1255,11 @@
 	        for (unsigned int qp=0; qp<n_qp; qp++)
 	          {
 	            // solution at the quadrature point
-	            Number fineval = fptr(xyz_values[qp],
-                                          parameters,
-                                          system.name(),
-                                          system.variable_name(var));
+	            Number fineval = (*f)(var_component,xyz_values[qp]);
 	            // solution grad at the quadrature point
 	            Gradient finegrad;
                     if (cont == C_ONE)
-                      finegrad = gptr(xyz_values[qp], parameters,
-                                      system.name(),
-				      system.variable_name(var));
+                      finegrad = (*g)(var_component,xyz_values[qp]);
 
                     // Form side projection matrix
                     for (unsigned int sidei=0, freei=0;
@@ -1354,16 +1341,11 @@
 	  for (unsigned int qp=0; qp<n_qp; qp++)
 	    {
 	      // solution at the quadrature point
-	      Number fineval = fptr(xyz_values[qp],
-                                    parameters,
-                                    system.name(),
-                                    system.variable_name(var));
+	      Number fineval = (*f)(var_component,xyz_values[qp]);
 	      // solution grad at the quadrature point
 	      Gradient finegrad;
               if (cont == C_ONE)
-                finegrad = gptr(xyz_values[qp], parameters,
-                                system.name(),
-				system.variable_name(var));
+                finegrad = (*g)(var_component,xyz_values[qp]);
 
               // Form interior projection matrix
               for (unsigned int i=0, freei=0; i != n_dofs; ++i)
Index: src/systems/system.C
===================================================================
--- src/systems/system.C	(revision 5063)
+++ src/systems/system.C	(working copy)
@@ -81,6 +81,7 @@
   _mesh                             (es.get_mesh()),
   _sys_name                         (name),
   _sys_number                       (number),
+  _variable_first_scalar            (1, 0),
   _active                           (true),
   _solution_projection              (true),
   _basic_system_only                (false),
@@ -182,6 +183,9 @@
 
   _variable_numbers.clear();
 
+  _variable_first_scalar.clear();
+  _variable_first_scalar.push_back(0);
+
   _dof_map->clear ();
 
   solution->clear ();
@@ -1031,6 +1035,12 @@
 
   _variable_numbers[var] = curr_n_vars;
 
+
+  libmesh_assert(_variable_first_scalar.size() == curr_n_vars+1);
+  _variable_first_scalar.push_back
+    (_variable_first_scalar[curr_n_vars] +
+     variable(curr_n_vars).n_components());
+
   // Add the variable to the _dof_map
   _dof_map->add_variable (_variables.back());