  #ifdef DEBUG
    {
      // Make sure we didn't miss any nodes
-@@ -1713,6 +1724,8 @@
+     MeshTools::libmesh_assert_valid_procids<Node>(mesh);
+-
++    
+     MeshBase::node_iterator       node_it  = mesh.local_nodes_begin();
+     const MeshBase::node_iterator node_end = mesh.local_nodes_end();
+     for (; node_it != node_end; ++node_it)
+       {
+ 	Node *obj = *node_it;
+ 	libmesh_assert(obj);
+-	unsigned int n_variables = obj->n_vars(this->sys_number());
+-	for (unsigned int v=0; v != n_variables; ++v)
++	unsigned int n_var_g = obj->n_var_groups(this->sys_number());
++	for (unsigned int vg=0; vg != n_var_g; ++vg)
+ 	  {
+-	    unsigned int n_comp =
+-	      obj->n_comp(this->sys_number(), v);
+-	    unsigned int first_dof = n_comp ?
+-	      obj->dof_number(this->sys_number(), v, 0) : 0;
++	    unsigned int n_comp_g =
++	      obj->n_comp_group(this->sys_number(), vg);
++	    unsigned int first_dof = n_comp_g ?
++	      obj->vg_dof_base(this->sys_number(), vg) : 0;
+ 	    libmesh_assert_not_equal_to (first_dof, DofObject::invalid_id);
+ 	  }
+       }
+@@ -1717,6 +1696,8 @@
    STOP_LOG("dof_indices()", "DofMap");
  }
  

   libmesh_assert_less (i, _variables.size());
Index: include/base/dof_object.h
===================================================================
--- include/base/dof_object.h	(revision 6500)
+++ include/base/dof_object.h	(working copy)
@@ -67,6 +67,7 @@
 

 
     for ( ; elem_it != elem_end; ++elem_it)
-      (*elem_it)->set_n_vars(this->sys_number(),n_var);
+      (*elem_it)->set_n_vars_per_group(sys_num, n_vars_per_group);
   }
 
 
@@ -526,18 +517,20 @@
 
   //------------------------------------------------------------
   // Next allocate space for the DOF indices
-  for (unsigned int var=0; var<this->n_variables(); var++)
+  for (unsigned int vg=0; vg<n_var_groups; vg++)
     {
-      const Variable &var_description =	this->variable(var);
-      const FEType& base_fe_type              = this->variable_type(var);
+      const VariableGroup &vg_description = this->variable_group(vg);
+      
+      const unsigned int n_var_in_group = vg_description.n_variables();
+      const FEType& base_fe_type        = vg_description.type();
 
       // Don't need to loop over elements for a SCALAR variable
       // Just increment _n_SCALAR_dofs
       if(base_fe_type.family == SCALAR)
-      {
-        this->_n_SCALAR_dofs += base_fe_type.order;
-        continue;
-      }
+	{
+	  this->_n_SCALAR_dofs += base_fe_type.order*n_var_in_group;
+	  continue;
+	}
 
       // This should be constant even on p-refined elements
       const bool extra_hanging_dofs =
@@ -554,13 +547,13 @@
 	  libmesh_assert(elem);
 
 	  // Skip the numbering if this variable is
-	  // not active on this element's subdoman
-	  if (!var_description.active_on_subdomain(elem->subdomain_id()))
+	  // not active on this element's subdomain
+	  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
 	    continue;
 
 	  const ElemType type = elem->type();
           const unsigned int dim = elem->dim();
-
+	  
           FEType fe_type = base_fe_type;
 
 #ifdef LIBMESH_ENABLE_AMR
@@ -585,7 +578,7 @@
 
                   libmesh_error();
                 }
-
+	      
               libMesh::err
                 << "WARNING: Finite element "
                 << Utility::enum_to_string(base_fe_type.family)
@@ -611,7 +604,7 @@
 	      if (elem->is_vertex(n))
 	        {
 	          const unsigned int old_node_dofs =
-	            node->n_comp(this->sys_number(), var);
+	            node->n_comp_group(sys_num, vg);
 
 		  const unsigned int vertex_dofs =
 		    std::max(FEInterface::n_dofs_at_node(dim, fe_type,
@@ -621,12 +614,13 @@
 		  // Some discontinuous FEs have no vertex dofs
 		  if (vertex_dofs > old_node_dofs)
 		    {
-		      node->set_n_comp(this->sys_number(), var,
+		      node->set_n_comp(sys_num, vg,
 				       vertex_dofs);
+		      
 		      // Abusing dof_number to set a "this is a
 		      // vertex" flag
-		      node->set_dof_number(this->sys_number(),
-					   var, 0, vertex_dofs);
+		      node->set_vg_dof_base(sys_num, vg,
+					    vertex_dofs);
 		    }
 	        }
 	    }
@@ -641,8 +635,8 @@
 	  libmesh_assert(elem);
 
 	  // Skip the numbering if this variable is
-	  // not active on this element's subdoman
-	  if (!var_description.active_on_subdomain(elem->subdomain_id()))
+	  // not active on this element's subdomain
+	  if (!vg_description.active_on_subdomain(elem->subdomain_id()))
 	    continue;
 
 	  const ElemType type = elem->type();
@@ -658,10 +652,10 @@
 	      Node* node = elem->get_node(n);
 
 	      const unsigned int old_node_dofs =
-	        node->n_comp(this->sys_number(), var);
+	        node->n_comp_group(sys_num, vg);
 
               const unsigned int vertex_dofs = old_node_dofs?
-                node->dof_number (this->sys_number(),var,0):0;
+                node->vg_dof_base (sys_num,vg):0;
 
 	      const unsigned int new_node_dofs =
 		FEInterface::n_dofs_at_node(dim, fe_type, type, n);
@@ -670,7 +664,7 @@
 	      if (elem->is_vertex(n))
 	        {
 		  libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
-		  libmesh_assert_greater_equal (vertex_dofs, new_node_dofs);
+		  libmesh_assert_greater_equal (vertex_dofs,   new_node_dofs);
 		}
 	      // We need to allocate the rest
 	      else
@@ -679,13 +673,13 @@
 		  // dofs, so we just give it edge or face dofs
 		  if (!old_node_dofs)
                     {
-		      node->set_n_comp(this->sys_number(), var,
+		      node->set_n_comp(sys_num, vg,
 				       new_node_dofs);
 		      // Abusing dof_number to set a "this has no
 		      // vertex dofs" flag
                       if (new_node_dofs)
-		        node->set_dof_number(this->sys_number(),
-					     var, 0, 0);
+		        node->set_vg_dof_base(sys_num, vg,
+					      0);
                     }
 
 		  // If this has dofs, but has no vertex dofs,
@@ -695,10 +689,11 @@
 		    {
                       if (new_node_dofs > old_node_dofs)
                         {
-		          node->set_n_comp(this->sys_number(), var,
+		          node->set_n_comp(sys_num, vg,
 				           new_node_dofs);
-		          node->set_dof_number(this->sys_number(),
-					       var, 0, vertex_dofs);
+			  
+		          node->set_vg_dof_base(sys_num, vg,
+						vertex_dofs);
                         }
 		    }
 		  // If this is another element's vertex,
@@ -708,10 +703,11 @@
 		    {
                       if (new_node_dofs > old_node_dofs - vertex_dofs)
                         {
-		          node->set_n_comp(this->sys_number(), var,
+		          node->set_n_comp(sys_num, vg,
 				           vertex_dofs + new_node_dofs);
-		          node->set_dof_number(this->sys_number(),
-					       var, 0, vertex_dofs);
+			  
+		          node->set_vg_dof_base(sys_num, vg,
+						vertex_dofs);
                         }
 		    }
 		  // If this is another element's vertex, add any
@@ -721,10 +717,11 @@
 		      libmesh_assert_greater_equal (old_node_dofs, vertex_dofs);
                       if (new_node_dofs > old_node_dofs)
                         {
-		          node->set_n_comp(this->sys_number(), var,
+		          node->set_n_comp(sys_num, vg,
 				           new_node_dofs);
-		          node->set_dof_number(this->sys_number(),
-					       var, 0, vertex_dofs);
+			  
+		          node->set_vg_dof_base (sys_num, vg,
+						 vertex_dofs);
                         }
 		    }
 		}
@@ -734,32 +731,15 @@
 			  FEInterface::n_dofs_per_elem(dim, fe_type,
 						       type);
 
-	  elem->set_n_comp(this->sys_number(), var, dofs_per_elem);
+	  elem->set_n_comp(sys_num, vg, dofs_per_elem);
 
 	}
-    }
+    } // end loop over variable groups
 
   // Calling DofMap::reinit() by itself makes little sense,
   // so we won't bother with nonlocal DofObjects.
   // Those will be fixed by distribute_dofs
-/*
-  //------------------------------------------------------------
-  // At this point, all n_comp and dof_number values on local
-  // DofObjects should be correct, but a ParallelMesh might have
-  // incorrect values on non-local DofObjects.  Let's request the
-  // correct values from each other processor.
 
-  this->set_nonlocal_n_comps(mesh.nodes_begin(),
-                             mesh.nodes_end(),
-                             mesh,
-                             &DofMap::node_ptr);
-
-  this->set_nonlocal_n_comps(mesh.elements_begin(),
-                             mesh.elements_end(),
-                             mesh,
-                             &DofMap::elem_ptr);
-*/
-
   //------------------------------------------------------------
   // Finally, clear all the current DOF indices
   // (distribute_dofs expects them cleared!)
@@ -799,9 +779,9 @@
   //_dof_coupling->clear();
 
   _variables.clear();
