  
  
  
-+const VariableGroup & DofMap::variable_group (const unsigned int g) const
-+{
-+  libmesh_assert_less (g, _variable_groups.size());
-+
-+  return _variable_groups[g];
-+}
-+
-+
-+
- const Variable & DofMap::variable (const unsigned int c) const
+-Order DofMap::variable_order (const unsigned int c) const
+-{
+-  libmesh_assert_less (c, _variables.size());
+-
+-  return _variables[c].type().order;
+-}
+-
+-
+-
+-const FEType& DofMap::variable_type (const unsigned int c) const
+-{
+-  libmesh_assert_less (c, _variables.size());
+-
+-  return _variables[c].type();
+-}
+-
+-
+-
+ void DofMap::attach_matrix (SparseMatrix<Number>& matrix)
  {
-   libmesh_assert_less (c, _variables.size());
-@@ -436,10 +456,15 @@
+   parallel_only();
+@@ -436,10 +420,17 @@
  
    START_LOG("reinit()", "DofMap");
  
 -  //this->clear();
 +  const unsigned int
-+    n_vars_tot   = this->n_variables(),
++    sys_num      = this->sys_number(),
 +    n_var_groups = this->n_variable_groups();
  
 -  const unsigned int n_var = this->n_variables();
 -
++  // The DofObjects need to know how many variable groups we have, and
++  // how many variables there are in each group.
 +  std::vector<unsigned int> n_vars_per_group; /**/ n_vars_per_group.reserve (n_var_groups);
 +  
 +  for (unsigned int vg=0; vg<n_var_groups; vg++)

 
     for ( ; elem_it != elem_end; ++elem_it)
-      (*elem_it)->set_n_vars(this->sys_number(),n_var);
+      (*elem_it)->set_n_vars_per_group(this->sys_number(), n_vars_per_group);
   }
 
 
@@ -799,9 +824,9 @@
   //_dof_coupling->clear();
 
   _variables.clear();

+  const unsigned int
+    start_idx_sys = this->start_idx(s),
+    nvg           = _idx_buf[start_idx_sys] / ncv_magic;
 
+
+  return nvg;
+}
+
+
+
 inline
+unsigned int DofObject::n_vars(const unsigned int s) const
+{
+  libmesh_assert_less (s, this->n_systems());
+
+  const unsigned int nvg = this->n_var_groups(s);
+    
+  unsigned int val=0;
+
+  for (unsigned int vg=0; vg<nvg; vg++)
+    val += this->n_vars(s,vg);
+
+  return val;
+}
+
+
+
+
+inline
 unsigned int DofObject::n_comp(const unsigned int s,
 			       const unsigned int var) const
 {
-  libmesh_assert_less (s, this->n_systems());
+  libmesh_assert_less (s,   this->n_systems());
   libmesh_assert_less (var, this->n_vars(s));
 
-# ifdef DEBUG
-  // Does this ever happen?  I doubt it... 3/7/2003 (BSK)
-  if (var >= this->n_vars(s))
-    {
-      libMesh::err << "s=" << s << ", var=" << var << std::endl
-		    << "this->n_vars(s)=" << this->n_vars(s) << std::endl
-		    << "this->n_systems()=" << this->n_systems() << std::endl;
-      libmesh_error();
-    }
-# endif
+  return this->n_comp_group(s,this->var_to_vg(s,var));
+}
 
+
+
+
+inline
+unsigned int DofObject::n_comp_group(const unsigned int s,
+				     const unsigned int vg) const
+{
+  libmesh_assert_less (s,  this->n_systems());
+  libmesh_assert_less (vg, this->n_var_groups(s));
+
   const unsigned int
     start_idx_sys = this->start_idx(s);
 
-  libmesh_assert_less ((start_idx_sys + 2*var), _idx_buf.size());
+  libmesh_assert_less ((start_idx_sys + 2*vg), _idx_buf.size());
 
-  return _idx_buf[start_idx_sys + 2*var];
+  return (_idx_buf[start_idx_sys + 2*vg] % ncv_magic);
 }
 
 
@@ -570,17 +690,20 @@
 				   const unsigned int var,
 				   const unsigned int comp) const
 {
-  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());
 
   const unsigned int
-    base_idx = _idx_buf[start_idx_sys + 2*var + 1];
+    base_idx = _idx_buf[start_idx_sys + 2*vg + 1];
 
   // if the first component is invalid, they
   // are all invalid
@@ -590,7 +713,7 @@
   // otherwise the index is the first component
   // index augemented by the component number
   else
-    return (base_idx + comp);
+    return (base_idx + vig*ncg + comp);
 }
 
 
@@ -639,6 +762,80 @@
 }
 
 
+
+inline
+void DofObject::set_vg_dof_base(const unsigned int s,
+				const unsigned int vg,
+				const unsigned int db)
+{
+  libmesh_assert_less (s,    this->n_systems());
+  libmesh_assert_less (vg,   this->n_var_groups(s));
+
+  const unsigned int
+    start_idx_sys = this->start_idx(s);
+
+  libmesh_assert_less ((start_idx_sys + 2*vg + 1), _idx_buf.size());
+
+  _idx_buf[start_idx_sys + 2*vg + 1] = db;
+
+  libmesh_assert_equal_to (this->vg_dof_base(s,vg), db);
+}
+
+
+
+inline
+unsigned int DofObject::vg_dof_base(const unsigned int s,
+				    const unsigned int vg) const
+{
+  libmesh_assert_less (s,    this->n_systems());
+  libmesh_assert_less (vg,   this->n_var_groups(s));
+
+  const unsigned int
+    start_idx_sys = this->start_idx(s);
+
+  libmesh_assert_less ((start_idx_sys + 2*vg + 1), _idx_buf.size());
+
+  return _idx_buf[start_idx_sys + 2*vg + 1];
+}
+
+
+
+inline
+unsigned int DofObject::var_to_vg (const unsigned int s,
