libMesh: A C++ Finite Element Library

You are proposing 

> return _var_first_local_df[2*var]

and

> return _var_first_local_df[2*var+1]

??

That is not consistent with the way the bounds are stored.  As I mentioned,
consider the three-variable system with unknowns (u,v,w).  For simplicity
say there are 150 dofs total, split 50-50-50.  The the _var_first_local_df
array would look like this:

_var_first_local_df[] = {0, 50, 100, 150};


The the u variable (0) is in the range

[_var_first_local_df[0], _var_first_local_df[1]) = [0,50)

the the v variable (1) is in the range

[_var_first_local_df[1], _var_first_local_df[2]) = [50,100)

and the the u variable (2) is in the range

[_var_first_local_df[2], _var_first_local_df[3]) = [100,150)

What you propose would assign the upper bound of w to
_var_first_local_df[2*2+1], which is out-of-bounds.

-Ben