 +   * ncv = ncv_magic*nv + nc
 +   * where nv are the number of identical variables of a given type,
 +   * and nc is the number of components for this set of variables.
-+   */
-+  static const index_t ncv_magic = 1000;
++   *
++   * It is hoped that by setting this to a power of two, an optimizing compiler
++   * will recgnize later that  #/ncv_magic is simply a bitshift
++   */
++  static const index_t ncv_magic = 256; // = 2^8, in case we want to manually bitshift
 +  
 +  /**
     * The starting index for system \p s.
     */
    unsigned int start_idx(const unsigned int s) const;
-@@ -527,6 +565,33 @@
+@@ -389,15 +473,15 @@
+   if (sys_num >= this->n_systems())
+     {
+       for (unsigned int s=0; s<this->n_systems(); s++)
+-        for (unsigned int v=0; v<this->n_vars(s); v++)
+-	  if (this->n_comp(s,v))
+-	    this->set_dof_number(s,v,0,invalid_id);
++        for (unsigned int vg=0; vg<this->n_var_groups(s); vg++)
++	  if (this->n_comp_group(s,vg))
++	    this->set_vg_dof_base(s,vg,invalid_id);
+     }
+   // ...otherwise invalidate the dofs for all systems
+   else
+-    for (unsigned int v=0; v<this->n_vars(sys_num); v++)
+-      if (this->n_comp(sys_num,v))
+-        this->set_dof_number(sys_num,v,0,invalid_id);
++    for (unsigned int vg=0; vg<this->n_var_groups(sys_num); vg++)
++      if (this->n_comp_group(sys_num,vg))
++        this->set_vg_dof_base(sys_num,vg,invalid_id);
+ }
+ 
+ 
+@@ -527,7 +611,7 @@
  
  
  inline
+-unsigned int DofObject::n_vars(const unsigned int s) const
 +unsigned int DofObject::n_var_groups(const unsigned int s) const
-+{
-+  libmesh_assert_less (s, this->n_systems());
-+
-+  return (this->end_idx(s) - this->start_idx(s)) / 2;
-+}
-+
-+
-+
+ {
+   libmesh_assert_less (s, this->n_systems());
+ 
+@@ -536,31 +620,67 @@
+ 
+ 
+ 
 +inline
 +unsigned int DofObject::n_vars(const unsigned int s,
 +			       const unsigned int vg) const

 #ifdef LIBMESH_ENABLE_AMR
 
   //------------------------------------------------------------
@@ -510,14 +535,14 @@
     const MeshBase::node_iterator node_end = mesh.nodes_end();
 
     for ( ; node_it != node_end; ++node_it)
-      (*node_it)->set_n_vars(this->sys_number(),n_var);
+      (*node_it)->set_n_vars_per_group(this->sys_number(), n_vars_per_group);
 
     // All the elements
     MeshBase::element_iterator       elem_it  = mesh.elements_begin();

   // 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?
