     * An invaild \p id to distinguish an uninitialized \p DofObject
     */
    static const unsigned int invalid_id = libMesh::invalid_uint;
-@@ -298,8 +336,22 @@
- 
+@@ -296,10 +342,29 @@
+    */
+   void pack_indexing(std::back_insert_iterator<std::vector<int> > target) const;
+ 
++  /**
++   * Print our buffer for debugging.
++   */
++  void debug_buffer () const;
++
  private:
  
 +  /**

 } // namespace libMesh
 
 #endif // #define __variable_h__
Index: src/systems/system.C
===================================================================
--- src/systems/system.C	(revision 6517)
+++ 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);
 }
 
 
Index: src/base/dof_object.C
===================================================================
--- src/base/dof_object.C	(revision 6517)
+++ 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(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,18 +293,23 @@
     _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
+  for (unsigned int vg=0; vg<this->n_var_groups(s); vg++)
+    libmesh_assert_equal_to (this->n_vars(s,vg), nvpg[vg]);
+   
   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)
@@ -302,29 +326,42 @@
 
 
 void DofObject::set_n_comp(const unsigned int s,
-			   const unsigned int var,
+			   const unsigned int vg,
 			   const unsigned int ncomp)
 {
-  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 (vg, this->n_var_groups(s));
+  
   // Check for trivial return
-  if (ncomp == this->n_comp(s,var)) return;
+  if (ncomp == this->n_comp(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;
+  // 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 +371,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;
 
@@ -364,6 +406,7 @@
   libmesh_assert_equal_to (this->dof_number(s, var, comp), dn);
 }
 
+  
 
 unsigned int DofObject::packed_indexing_size() const
 {
@@ -377,6 +420,7 @@
 }
 
 
+  
 unsigned int DofObject::unpackable_indexing_size(std::vector<int>::const_iterator begin)
 {
 #ifdef LIBMESH_ENABLE_AMR
Index: src/base/dof_map.C
===================================================================
--- src/base/dof_map.C	(revision 6517)
+++ src/base/dof_map.C	(working copy)
@@ -128,11 +128,11 @@
 DofMap::DofMap(const unsigned int number) :

 } // namespace libMesh
 
 
Index: include/base/variable.h
===================================================================
--- include/base/variable.h	(revision 6518)
+++ include/base/variable.h	(working copy)
@@ -19,12 +19,14 @@
 #define __variable_h__
 
 // Local Includes
+#include "libmesh/libmesh_common.h"
 #include "libmesh/fe_type.h"
 #include "libmesh/id_types.h"
 
 // C++ includes
 #include <set>
 #include <string>
+#include <vector>
 
 namespace libMesh {
 
@@ -109,7 +111,7 @@
    * empty \p _active_subdomains container as active everywhere, i.e.
    * for all subdomains.
    */
