  void DofMap::SCALAR_dof_indices (std::vector<unsigned int>& di,
  			         const unsigned int vn,
  #ifdef LIBMESH_ENABLE_AMR
+Index: src/base/dof_object.C
+===================================================================
+--- src/base/dof_object.C	(revision 6518)
++++ src/base/dof_object.C	(working copy)
+@@ -55,8 +55,12 @@
+ 
+   for (unsigned int s=0; s<this->n_systems(); s++)
+     {
+-      libmesh_assert_equal_to (this->n_vars(s), dof_obj.n_vars(s));
++      libmesh_assert_equal_to (this->n_vars(s),       dof_obj.n_vars(s));
++      libmesh_assert_equal_to (this->n_var_groups(s), dof_obj.n_var_groups(s));
+ 
++      for (unsigned int vg=0; vg<this->n_var_groups(s); vg++)
++	libmesh_assert_equal_to (this->n_vars(s,vg), dof_obj.n_vars(s,vg));
++      
+       for (unsigned int v=0; v<this->n_vars(s); v++)
+ 	{
+ 	  libmesh_assert_equal_to (this->n_comp(s,v), dof_obj.n_comp(s,v));
+@@ -91,8 +95,12 @@
+ 
+   for (unsigned int s=0; s<this->n_systems(); s++)
+     {
+-      libmesh_assert_equal_to (this->n_vars(s), dof_obj.n_vars(s));
++      libmesh_assert_equal_to (this->n_vars(s),       dof_obj.n_vars(s));
++      libmesh_assert_equal_to (this->n_var_groups(s), dof_obj.n_var_groups(s));
+ 
++      for (unsigned int vg=0; vg<this->n_var_groups(s); vg++)
++	libmesh_assert_equal_to (this->n_vars(s,vg), dof_obj.n_vars(s,vg));
++      
+       for (unsigned int v=0; v<this->n_vars(s); v++)
+ 	{
+ 	  libmesh_assert_equal_to (this->n_comp(s,v), dof_obj.n_comp(s,v));
+@@ -161,8 +169,11 @@
+   // check that all systems now exist and that they have 0 size
+   libmesh_assert_equal_to (ns, this->n_systems());
+   for (unsigned int s=0; s<this->n_systems(); s++)
+-    libmesh_assert_equal_to (this->n_vars(s), 0);
+-
++    {
++      libmesh_assert_equal_to (this->n_vars(s),       0);
++      libmesh_assert_equal_to (this->n_var_groups(s), 0);
++    }
++  
+ #endif
+ }
+ 
+@@ -195,15 +206,20 @@
+ 
+   libmesh_assert_equal_to (this->n_systems(), (ns_orig+1));
+   libmesh_assert_equal_to (this->n_vars(ns_orig), 0);
++  libmesh_assert_equal_to (this->n_var_groups(ns_orig), 0);
+ }
+ 
+ 
+ 
+-void DofObject::set_n_vars(const unsigned int s,
+-			   const unsigned int nvars)
++void DofObject::set_n_vars_per_group(const unsigned int s,
++				     const std::vector<unsigned int> &nvpg)
+ {
++  
+   libmesh_assert_less (s, this->n_systems());
+ 
++  // number of varaible groups for this system - inferred
++  const unsigned int nvg = nvpg.size();
++  
+   // BSK - note that for compatibility with the previous implementation
+   // calling this method when (nvars == this->n_vars()) requires that
+   // we invalidate the DOF indices and set the number of components to 0.
+@@ -212,10 +228,13 @@
+   // even in the case of (nvars == this->n_vars()), resulting in n_comp(s,v)
+   // implicitly becoming 0 regardless of any previous value.
+   // quick return?
+-  if (nvars == this->n_vars(s))
++  if (nvg == this->n_var_groups(s))
+     {
+-      for (unsigned int v=0; v<nvars; v++)
+-	this->set_n_comp(s,v,0);
++      for (unsigned int vg=0; vg<nvg; vg++)
++	{
++	  this->set_n_comp_group(s,vg,0);
++	  libmesh_assert_equal_to (this->n_vars(s,vg), nvpg[vg]);
++	}
+       return;
+     }
+ 
+@@ -226,13 +245,13 @@
+   old_system_sizes.reserve(this->n_systems());
+ 
+   for (unsigned int s_ctr=0; s_ctr<this->n_systems(); s_ctr++)
+-    old_system_sizes.push_back(this->n_vars(s_ctr));
++    old_system_sizes.push_back(this->n_var_groups(s_ctr));
+ #endif
+ 
+   // remove current indices if we have some
+-  if (this->n_vars(s) != 0)
++  if (this->n_var_groups(s) != 0)
+     {
+-      const unsigned int old_nvars_s = this->n_vars(s);
++      const unsigned int old_nvg_s = this->n_var_groups(s);
+ 
+       DofObject::index_buffer_t::iterator
+ 	it  = _idx_buf.begin(),
+@@ -243,30 +262,30 @@
+       _idx_buf.erase(it,end);
+ 
+       for (unsigned int ctr=(s+1); ctr<this->n_systems(); ctr++)
+-	_idx_buf[ctr] -= 2*old_nvars_s;
++	_idx_buf[ctr] -= 2*old_nvg_s;
+     }
+ 
+   // better not have any now!
+-  libmesh_assert_equal_to (this->n_vars(s), 0);
++  libmesh_assert_equal_to (this->n_var_groups(s), 0);
+ 
+   // had better not screwed up any of our sizes!
+ #ifdef DEBUG
+   for (unsigned int s_ctr=0; s_ctr<this->n_systems(); s_ctr++)
+     if (s_ctr != s)
+-      libmesh_assert_equal_to (this->n_vars(s_ctr), old_system_sizes[s_ctr]);
++      libmesh_assert_equal_to (this->n_var_groups(s_ctr), old_system_sizes[s_ctr]);
+ #endif
+ 
+   // OK, if the user requested 0 that is what we have
+-  if (nvars == 0)
++  if (nvg == 0)
+     return;
+ 
+   {
+     // array to hold new indices
+-    DofObject::index_buffer_t var_idxs(2*nvars);
+-    for (unsigned int v=0; v<nvars; v++)
+-      {
+-	var_idxs[2*v    ] = 0;
+-	var_idxs[2*v + 1] = invalid_id - 1;
++    DofObject::index_buffer_t var_idxs(2*nvg);
++    for (unsigned int vg=0; vg<nvg; vg++)
++      {	
++	var_idxs[2*vg    ] = ncv_magic*nvpg[vg] + 0;
++	var_idxs[2*vg + 1] = invalid_id - 1;
+       }
+ 
+     DofObject::index_buffer_t::iterator it = _idx_buf.begin();
+@@ -274,28 +293,38 @@
+     _idx_buf.insert(it, var_idxs.begin(), var_idxs.end());
+ 
+     for (unsigned int ctr=(s+1); ctr<this->n_systems(); ctr++)
+-      _idx_buf[ctr] += 2*nvars;
++      _idx_buf[ctr] += 2*nvg;
+ 
+     // resize _idx_buf to fit so no memory is wasted.
+     DofObject::index_buffer_t(_idx_buf).swap(_idx_buf);
+   }
+ 
+   // that better had worked.  Assert stuff.
+-  libmesh_assert_equal_to (nvars, this->n_vars(s));
++  libmesh_assert_equal_to (nvg, this->n_var_groups(s));
++  
++#ifdef DEBUG
+ 
+-#ifdef DEBUG
++  // std::cout << " [ ";
++  // for (unsigned int i=0; i<_idx_buf.size(); i++)
++  //   std::cout << _idx_buf[i] << " ";
++  // std::cout << "]\n";
++
++  libmesh_assert_equal_to (this->n_var_groups(s), nvpg.size());
++
++  for (unsigned int vg=0; vg<this->n_var_groups(s); vg++)
++    {
++      libmesh_assert_equal_to (this->n_vars(s,vg), nvpg[vg]);
++      libmesh_assert_equal_to (this->n_comp_group(s,vg), 0);
++    }
++
+   for (unsigned int v=0; v<this->n_vars(s); v++)
+     libmesh_assert_equal_to (this->n_comp(s,v), 0);
++  
+   // again, all other system sizes shoudl be unchanged!
+   for (unsigned int s_ctr=0; s_ctr<this->n_systems(); s_ctr++)
+     if (s_ctr != s)
+-      libmesh_assert_equal_to (this->n_vars(s_ctr), old_system_sizes[s_ctr]);
++      libmesh_assert_equal_to (this->n_var_groups(s_ctr), old_system_sizes[s_ctr]);
+ 
+-//   std::cout << " [ ";
+-//   for (unsigned int i=0; i<_idx_buf.size(); i++)
+-//     std::cout << _idx_buf[i] << " ";
+-//   std::cout << "]\n";
+-
+ #endif
+ }
+ 
+@@ -305,26 +334,57 @@
+ 			   const unsigned int var,
+ 			   const unsigned int ncomp)
+ {
+-  libmesh_assert_less (s, this->n_systems());
++  libmesh_assert_less (s,   this->n_systems());
+   libmesh_assert_less (var, this->n_vars(s));
+ 
++  this->set_n_comp_group(s, this->var_to_vg(s,var), ncomp);
++}
++
++
++
++void DofObject::set_n_comp_group(const unsigned int s,
++				 const unsigned int vg,
++				 const unsigned int ncomp)
++{
++  libmesh_assert_less (s,  this->n_systems());
++  libmesh_assert_less (vg, this->n_var_groups(s));
++  
+   // Check for trivial return
+-  if (ncomp == this->n_comp(s,var)) return;
++  if (ncomp == this->n_comp_group(s,vg)) return;
+ 
++#ifndef NDEBUG
++  if (ncomp >= ncv_magic)
++    {
++      const index_t ncvm = ncv_magic;
++      libMesh::err << "ERROR: ncomp must be less than DofObject::ncv_magic!\n"
++		   << "ncomp = " << ncomp << ", ncv_magic = " << ncvm
++		   << "\nrecompile and try again!\n";
++      libmesh_error();
++    }
++#endif
++  
+   const unsigned int
+     start_idx_sys = this->start_idx(s),
+-    base_offset  = start_idx_sys + 2*var;
++    n_vars_group  = this->n_vars(s,vg),
++    base_offset   = start_idx_sys + 2*vg;
+ 
+   libmesh_assert_less ((base_offset + 1), _idx_buf.size());
+ 
+-  // set the number of components
+-  _idx_buf[base_offset] = ncomp;
++  // if (ncomp)
++  //   std::cout << "s,vg,ncomp=" 
++  // 	      << s  << ","
++  // 	      << vg << ","
++  // 	      << ncomp << '\n';
+ 
++  // set the number of components, maintaining the number
++  // of variables in the group
++  _idx_buf[base_offset] = ncv_magic*n_vars_group + ncomp;
++
+   // We use (invalid_id - 1) to signify no
+   // components for this object
+   _idx_buf[base_offset + 1] = (ncomp == 0) ? invalid_id - 1 : invalid_id;
+ 
+-  libmesh_assert_equal_to (ncomp, this->n_comp(s,var));
++  libmesh_assert_equal_to (ncomp, this->n_comp(s,vg));
+ }
+ 
+ 
+@@ -334,22 +394,27 @@
+ 			       const unsigned int comp,
+ 			       const unsigned int dn)
+ {
+-  libmesh_assert_less (s, this->n_systems());
+-  libmesh_assert_less (var, this->n_vars(s));
++  libmesh_assert_less (s,    this->n_systems());
++  libmesh_assert_less (var,  this->n_vars(s));
+   libmesh_assert_less (comp, this->n_comp(s,var));
+ 
+   const unsigned int
++    vg            = this->var_to_vg(s,var),
++    ncg           = this->n_comp_group(s,vg),
++    vig           = this->system_var_to_vg_var(s,vg,var),
+     start_idx_sys = this->start_idx(s);
+ 
+-  libmesh_assert_less ((start_idx_sys + 2*var + 1), _idx_buf.size());
++  libmesh_assert_less ((start_idx_sys + 2*vg + 1), _idx_buf.size());
+ 
+   unsigned int
+-    &base_idx(_idx_buf[start_idx_sys + 2*var + 1]);
++    &base_idx(_idx_buf[start_idx_sys + 2*vg + 1]);
+ 
+-  //We intend to change all dof numbers together or not at all
+-  if (comp)
++  // We intend to change all dof numbers together or not at all
++  if (comp || vig)
+     libmesh_assert ((dn == invalid_id && base_idx == invalid_id) ||
+-		    (dn == base_idx + comp));
++		    (dn == base_idx + vig*ncg + comp));
++  
++  // only explicitly store the base index for vig==0, comp==0
+   else
+     base_idx = dn;
+ 
+@@ -358,12 +423,12 @@
+ //   for (unsigned int i=0; i<_idx_buf.size(); i++)
+ //     std::cout << _idx_buf[i] << " ";
+ //   std::cout << "]\n";
+-
+ // #endif
+ 
+   libmesh_assert_equal_to (this->dof_number(s, var, comp), dn);
+ }
+ 
++  
+ 
+ unsigned int DofObject::packed_indexing_size() const
+ {
+@@ -377,6 +442,7 @@
+ }
+ 
+ 
++  
+ unsigned int DofObject::unpackable_indexing_size(std::vector<int>::const_iterator begin)
+ {
+ #ifdef LIBMESH_ENABLE_AMR
+@@ -454,4 +520,14 @@
+ 
+ 
+ 
++void DofObject::debug_buffer () const
++{
++  std::cout << " [ ";
++  for (unsigned int i=0; i<_idx_buf.size(); i++)
++    std::cout << _idx_buf[i] << " ";
++  std::cout << "]\n";
++}
++
++
++
+ } // namespace libMesh
+Index: src/systems/system.C
+===================================================================
+--- src/systems/system.C	(revision 6518)
++++ src/systems/system.C	(working copy)
+@@ -1073,53 +1073,86 @@
+ 			           const FEType& type,
+ 				   const std::set<subdomain_id_type> * const active_subdomains)
+ {
++  return this->add_variables (std::vector<std::string>(1, var),
++			      type,
++			      active_subdomains);
++}
++
++
++
++unsigned int System::add_variable (const std::string& var,
++			           const Order order,
++			           const FEFamily family,
++				   const std::set<subdomain_id_type> * const active_subdomains)
++{
++  return this->add_variables(std::vector<std::string>(1, var),
++			     FEType(order, family),
++			     active_subdomains);
++}
++
++
++
++unsigned int System::add_variables (const std::vector<std::string> &vars,
++				    const FEType& type,
++				    const std::set<subdomain_id_type> * const active_subdomains)
++{
++  const std::string var = vars[0];
++  
+   // Make sure the variable isn't there already
+   // or if it is, that it's the type we want
+-  for (unsigned int v=0; v<this->n_vars(); v++)
+-    if (this->variable_name(v) == var)
+-      {
+-	if (this->variable_type(v) == type)
+-	  return _variables[v].number();
++  for (unsigned int ov=0; ov<vars.size(); ov++)
++    for (unsigned int v=0; v<this->n_vars(); v++)
++      if (this->variable_name(v) == vars[ov])
++	{
++	  if (this->variable_type(v) == type)
++	    return _variables[v].number();
++	  
++	  libMesh::err << "ERROR: incompatible variable "
++		       << vars[ov]
++		       << " has already been added for this system!"
++		       << std::endl;
++	  libmesh_error();
++	}
+ 
+-	libMesh::err << "ERROR: incompatible variable "
+-		      << var
+-		      << " has already been added for this system!"
+-		      << std::endl;
+-	libmesh_error();
+-      }
+-
+   const unsigned int curr_n_vars = this->n_vars();
+ 
+   const unsigned int next_first_component = this->n_components();
+ 
+-  // Add the variable to the list
+-  _variables.push_back((active_subdomains == NULL) ?
+-		       Variable(var, curr_n_vars,
+-				next_first_component, type) :
+-		       Variable(var, curr_n_vars,
+-				next_first_component, type, *active_subdomains));
++  // Add the variable group to the list
++  _variable_groups.push_back((active_subdomains == NULL) ?
++			     VariableGroup(vars, curr_n_vars,
++					   next_first_component, type) :
++			     VariableGroup(vars, curr_n_vars,
++					   next_first_component, type, *active_subdomains));
+ 
+-  libmesh_assert_equal_to ((curr_n_vars+1), this->n_vars());
++  const VariableGroup &vg (_variable_groups.back());
++  
++  // Add each component of the group individually
++  for (unsigned int v=0; v<vars.size(); v++)
++    {
++      _variables.push_back (vg(v));
++      _variable_numbers[vars[v]] = curr_n_vars+v;  
++    }
++  
++  libmesh_assert_equal_to ((curr_n_vars+vars.size()), this->n_vars());
++    
++  // Add the variable group to the _dof_map
++  _dof_map->add_variable_group (vg);
+ 
+-  _variable_numbers[var] = curr_n_vars;
+-
+-  // Add the variable to the _dof_map
+-  _dof_map->add_variable (_variables.back());
+-
+   // Return the number of the new variable
+   return curr_n_vars;
+ }
+ 
+ 
+ 
+-unsigned int System::add_variable (const std::string& var,
+-			           const Order order,
+-			           const FEFamily family,
+-				   const std::set<subdomain_id_type> * const active_subdomains)
++unsigned int System::add_variables (const std::vector<std::string> &vars,
++				    const Order order,
++				    const FEFamily family,
++				    const std::set<subdomain_id_type> * const active_subdomains)
+ {
+-  return this->add_variable(var,
+-			    FEType(order, family),
+-			    active_subdomains);
++  return this->add_variables(vars,
++			     FEType(order, family),
++			     active_subdomains);
+ }
+ 
+ 
+@@ -1535,9 +1568,16 @@
+       << "    Type \""  << this->system_type() << "\"\n"
+       << "    Variables=";
+ 
+-  for (unsigned int vn=0; vn<this->n_vars(); vn++)
+-      out << "\"" << this->variable_name(vn) << "\" ";
++  for (unsigned int vg=0; vg<this->n_variable_groups(); vg++)
++    {
++      const VariableGroup &vg_description (this->variable_group(vg));
+ 
++      if (vg_description.n_variables() > 1) out << "{ ";
++      for (unsigned int vn=0; vn<vg_description.n_variables(); vn++)
++	out << "\"" << vg_description.name(vn) << "\" ";
++      if (vg_description.n_variables() > 1) out << "}";      
++    }
++
+   out << '\n';
+ 
+   out << "    Finite Element Types=";

 } // namespace libMesh
 
 
Index: include/base/dof_map.h
===================================================================
--- include/base/dof_map.h	(revision 6517)
+++ include/base/dof_map.h	(working copy)
@@ -309,13 +309,24 @@
     return *_n_oz;
   }
 
+  // /**
+  //  * Add an unknown of order \p order and finite element type
+  //  * \p type to the system of equations.
+  //  */
+  // void add_variable (const Variable &var);
+
   /**
-   * Add an unknown of order \p order and finite element type
+   * Add a group of unknowns of order \p order and finite element type
    * \p type to the system of equations.
    */
-  void add_variable (const Variable &var);
+  void add_variable_group (const VariableGroup &var_group);
 
   /**
+   * @returns the \p VariableGroup description object for group \p g.
+   */
+  const VariableGroup& variable_group (const unsigned int c) const;
+
+  /**
    * @returns the variable description object for variable \p c.
    */
   const Variable& variable (const unsigned int c) const;
@@ -326,23 +337,41 @@
   Order variable_order (const unsigned int c) const;
 
   /**
+   * @returns the approximation order for \p VariableGroup \p vg.
+   */
+  Order variable_group_order (const unsigned int vg) const;
+
+  /**
    * @returns the finite element type for variable \p c.
    */
   const FEType& variable_type (const unsigned int c) const;
 
   /**
+   * @returns the finite element type for \p VariableGroup \p vg.
+   */
+  const FEType& variable_group_type (const unsigned int vg) const;
+
+  /**
    * Returns the number of variables in the global solution vector. Defaults
    * to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier
    * Stokes (u,v,p), etc...
    */
+  unsigned int n_variable_groups() const
+  { return _variable_groups.size(); }
+
+  /**
+   * Returns the number of variables in the global solution vector. Defaults
+   * to 1, should be 1 for a scalar equation, 3 for 2D incompressible Navier
+   * Stokes (u,v,p), etc...
+   */
   unsigned int n_variables() const
   { return _variables.size(); }
 
   /**
    * @returns the total number of degrees of freedom in the problem.
    */
-
   unsigned int n_dofs() const { return _n_dfs; }
+  
   /**
    * @returns the number of SCALAR dofs.
    */
@@ -399,27 +428,6 @@
 #endif //LIBMESH_ENABLE_AMR
 
   /**
-   * Returns the first local degree of freedom index for variable \p var.
-   */
-  unsigned int variable_first_local_dof (const unsigned int var) const
-  {
-    libmesh_assert_less ((var+1), _var_first_local_df.size());
-    libmesh_assert_not_equal_to (_var_first_local_df[var], DofObject::invalid_id);
-    return _var_first_local_df[var];
-  }
-
-  /**
-   * Returns (one past) the last local degree of freedom index for variable \p var.
-   * Analogous to the end() member function of STL containers.
-   */
-  unsigned int variable_last_local_dof (const unsigned int var) const
-  {
-    libmesh_assert_less ((var+1), _var_first_local_df.size());
-    libmesh_assert_not_equal_to (_var_first_local_df[var+1], DofObject::invalid_id);
-    return _var_first_local_df[var+1];
-  }
-
-  /**
    * Fills the vector \p di with the global degree of freedom indices
    * for the element. If no variable number is specified then all
    * variables are returned.
@@ -979,14 +987,18 @@
   void add_constraints_to_send_list();
 
 #endif // LIBMESH_ENABLE_CONSTRAINTS
-
-
+  
   /**
    * The finite element type for each variable.
    */
   std::vector<Variable> _variables;
 
   /**
+   * The finite element type for each variable.
+   */
+  std::vector<VariableGroup> _variable_groups;
+
+  /**
    * The number of the system we manage DOFs for.
    */
   const unsigned int _sys_number;
@@ -1009,18 +1021,13 @@
   std::vector<unsigned int> _end_df;
 
   /**
-   * The first local DOF index for each variable in the \p System.
-   */
-  std::vector<unsigned int> _var_first_local_df;
-
-  /**
    * A list containing all the global DOF indicies that affect the
    * solution on my subdomain.
    */
   std::vector<unsigned int> _send_list;
 
   /**
-   *
+   * Funtion object to call to add extra entries to the sparsity pattern
    */
   AugmentSparsityPattern *_augment_sparsity_pattern;
 
@@ -1037,7 +1044,7 @@
   void * _extra_sparsity_context;
 
   /**
-   * 
+   * Function object to call to add extra entries to the send list
    */
   AugmentSendList *_augment_send_list;
 
@@ -1145,14 +1152,75 @@
 
 // ------------------------------------------------------------
 // Dof Map inline member functions
-
 inline
 unsigned int DofMap::sys_number() const
 {
   return _sys_number;
 }
 
+
+
 inline
+const VariableGroup & DofMap::variable_group (const unsigned int g) const
+{
+  libmesh_assert_less (g, _variable_groups.size());
+
+  return _variable_groups[g];
+}
+
+
+
+inline
+const Variable & DofMap::variable (const unsigned int c) const
+{
+  libmesh_assert_less (c, _variables.size());
+
+  return _variables[c];
+}
+
+
+
+inline
+Order DofMap::variable_order (const unsigned int c) const
+{
+  libmesh_assert_less (c, _variables.size());
+
+  return _variables[c].type().order;
+}
+
+
+
+inline
+Order DofMap::variable_group_order (const unsigned int vg) const
+{
+  libmesh_assert_less (vg, _variable_groups.size());
+
+  return _variable_groups[vg].type().order;
+}
+
+
+
+inline
+const FEType& DofMap::variable_type (const unsigned int c) const
+{
+  libmesh_assert_less (c, _variables.size());
+
+  return _variables[c].type();
+}
+
+
+
+inline
+const FEType& DofMap::variable_group_type (const unsigned int vg) const
+{
+  libmesh_assert_less (vg, _variable_groups.size());
+
+  return _variable_groups[vg].type();
+}
+
+
+
+inline
 bool DofMap::is_constrained_node (const Node*
 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
 node
Index: include/base/variable.h
===================================================================
--- include/base/variable.h	(revision 6517)
+++ include/base/variable.h	(working copy)
@@ -109,7 +109,7 @@
    * empty \p _active_subdomains container as active everywhere, i.e.
    * for all subdomains.
    */

         continue;
 
       MeshBase::element_iterator       elem_it  = mesh.active_local_elements_begin();
@@ -1115,7 +1112,7 @@
 
 	  // ... and only variables which are active on
 	  // on this element's subdomain
