diff --git a/framework/include/actions/AddArrayBCAction.h b/framework/include/actions/AddArrayBCAction.h
new file mode 100644
index 000000000000..0f3f13de2ca1
--- /dev/null
+++ b/framework/include/actions/AddArrayBCAction.h
@@ -0,0 +1,34 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#ifndef ADDARRAYBCACTION_H
+#define ADDARRAYBCACTION_H
+
+#include "MooseObjectAction.h"
+
+class AddArrayBCAction;
+
+template<>
+InputParameters validParams<AddArrayBCAction>();
+
+
+class AddArrayBCAction : public MooseObjectAction
+{
+public:
+  AddArrayBCAction(InputParameters params);
+
+  virtual void act() override;
+};
+
+#endif // ADDARRAYBCACTION_H
diff --git a/framework/include/base/ArrayMooseVariable.h b/framework/include/base/ArrayMooseVariable.h
index 7b90dadd9dea..d2d1428fa5db 100644
--- a/framework/include/base/ArrayMooseVariable.h
+++ b/framework/include/base/ArrayMooseVariable.h
@@ -62,7 +62,7 @@ class SystemBase;
 class ArrayMooseVariable : public MooseVariableBase
 {
 public:
-  ArrayMooseVariable(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count);
+  ArrayMooseVariable(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count);
   virtual ~ArrayMooseVariable();
 
   /**
@@ -205,10 +205,10 @@ class ArrayMooseVariable : public MooseVariableBase
 //   * Compute values at quadrature points for the neighbor
 //   */
 //  void computeNeighborValues();
-//  /**
-//   * Compute nodal values of this variable
-//   */
-//  void computeNodalValues();
+    /**
+     * Compute nodal values of this variable
+     */
+    void computeNodalValues();
 //  /**
 //   * Compute nodal values of this variable in the neighbor
 //   */
diff --git a/framework/include/base/Factory.h b/framework/include/base/Factory.h
index 41a728b78a85..1d3e790f780d 100644
--- a/framework/include/base/Factory.h
+++ b/framework/include/base/Factory.h
@@ -38,6 +38,7 @@
 #define registerKernel(name)                        registerObject(name)
 #define registerNodalKernel(name)                   registerObject(name)
 #define registerBoundaryCondition(name)             registerObject(name)
+#define registerArrayBoundaryCondition(name)        registerObject(name)
 #define registerAux(name)                           registerObject(name)
 #define registerAuxKernel(name)                     registerObject(name)
 #define registerMaterial(name)                      registerObject(name)
diff --git a/framework/include/base/MooseVariable.h b/framework/include/base/MooseVariable.h
index 435f766adff5..9960f301b9a9 100644
--- a/framework/include/base/MooseVariable.h
+++ b/framework/include/base/MooseVariable.h
@@ -41,7 +41,7 @@ template <typename T> class DenseVector;
 class MooseVariable : public MooseVariableBase
 {
 public:
-  MooseVariable(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind);
+  MooseVariable(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind);
   virtual ~MooseVariable();
 
   /**
diff --git a/framework/include/base/MooseVariableBase.h b/framework/include/base/MooseVariableBase.h
index 19a249fdccee..64b7233515d2 100644
--- a/framework/include/base/MooseVariableBase.h
+++ b/framework/include/base/MooseVariableBase.h
@@ -50,7 +50,7 @@ class SystemBase;
 class MooseVariableBase
 {
 public:
-  MooseVariableBase(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count = 1);
+  MooseVariableBase(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count = 1);
   virtual ~MooseVariableBase();
 
   /**
@@ -70,10 +70,16 @@ class MooseVariableBase
   SystemBase & sys() { return _sys; }
 
   /**
-   * Get the variable number
+   * Get the variable name
    */
   const std::string & name() const;
 
+  /**
+   * Get the names of all of the variables represented by this MoosevariableBase
+   * Note: this is usually built on the fly, so try not to call it too much!
+   */
+  const std::vector<std::string> names() const;
+
   /**
    * Kind of the variable (Nonlinear, Auxiliary, ...)
    */
@@ -134,9 +140,11 @@ class MooseVariableBase
   /**
    * The number of variables represented by this ArrayMooseVariable.
    */
-  unsigned int count() { return _count; }
+  unsigned int count() const { return _count; }
 
 protected:
+  /// The name of the variable
+  std::string _name;
   /// variable number (from libMesh)
   unsigned int _var_num;
   /// The FEType associated with this variable
diff --git a/framework/include/base/MooseVariableScalar.h b/framework/include/base/MooseVariableScalar.h
index 60e92289b4c9..713310fe41b1 100644
--- a/framework/include/base/MooseVariableScalar.h
+++ b/framework/include/base/MooseVariableScalar.h
@@ -30,7 +30,7 @@ template <typename T> class NumericVector;
 class MooseVariableScalar : public MooseVariableBase
 {
 public:
-  MooseVariableScalar(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind);
+  MooseVariableScalar(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind);
   virtual ~MooseVariableScalar();
 
   void reinit();
diff --git a/framework/include/base/NonlinearSystem.h b/framework/include/base/NonlinearSystem.h
index 9e2b6f5708bf..62057def66c9 100644
--- a/framework/include/base/NonlinearSystem.h
+++ b/framework/include/base/NonlinearSystem.h
@@ -33,7 +33,7 @@ class Predictor;
 class ElementDamper;
 class GeneralDamper;
 class IntegratedBC;
-class NodalBC;
+class NodalBCBase;
 class PresetNodalBC;
 class DGKernel;
 class InterfaceKernel;
@@ -529,7 +529,7 @@ class NonlinearSystem : public SystemTempl<TransientNonlinearImplicitSystem>,
   ///@{
   /// BoundaryCondition Warhouses
   MooseObjectWarehouse<IntegratedBC> _integrated_bcs;
-  MooseObjectWarehouse<NodalBC> _nodal_bcs;
+  MooseObjectWarehouse<NodalBCBase> _nodal_bcs;
   MooseObjectWarehouse<PresetNodalBC> _preset_nodal_bcs;
   ///@}
 
diff --git a/framework/include/base/SystemBase.h b/framework/include/base/SystemBase.h
index 1dbcd56e4636..8f1aa783b804 100644
--- a/framework/include/base/SystemBase.h
+++ b/framework/include/base/SystemBase.h
@@ -495,7 +495,7 @@ class SystemTempl : public SystemBase
     for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
     {
       //FIXME: we cannot refer fetype in libMesh at this point, so we will just make a copy in MooseVariableBase.
-      MooseVariable * var = new MooseVariable(var_num, type, *this, _subproblem.assembly(tid), _var_kind);
+      MooseVariable * var = new MooseVariable(var_name, var_num, type, *this, _subproblem.assembly(tid), _var_kind);
       var->scalingFactor(scale_factor);
       _vars[tid].add(var_name, var);
     }
@@ -513,8 +513,6 @@ class SystemTempl : public SystemBase
    */
   virtual void addArrayVariable(const std::string & var_name, const FEType & type, Real scale_factor, unsigned int count, const std::set< SubdomainID > * const active_subdomains = NULL)
   {
-    std::cout<<"addArrayVariable : "<<var_name<<std::endl;
-
     // Turn off automatic variable group identification so that we can be sure that this variable
     // group will be ordered exactly like it should be
     _sys.identify_variable_groups(false);
@@ -524,10 +522,8 @@ class SystemTempl : public SystemBase
     for (unsigned int i = 0; i < count; i++)
       var_names.push_back(var_name + "_" + std::to_string(i));
 
-    unsigned int var_num = _sys.add_variables(var_names, type, active_subdomains);
-
-    // FIXME: need to get this properly!
-    var_num = 0;
+    // The number returned by libMesh is the _last_ variable number... we want to hold onto the _first_
+    unsigned int var_num = _sys.add_variables(var_names, type, active_subdomains) - (count - 1);
 
     if (active_subdomains == NULL)
     {
@@ -542,7 +538,7 @@ class SystemTempl : public SystemBase
     for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
     {
       //FIXME: we cannot refer fetype in libMesh at this point, so we will just make a copy in MooseVariableBase.
-      ArrayMooseVariable * var = new ArrayMooseVariable(var_num, type, *this, _subproblem.assembly(tid), _var_kind, count);
+      ArrayMooseVariable * var = new ArrayMooseVariable(var_name, var_num, type, *this, _subproblem.assembly(tid), _var_kind, count);
       var->scalingFactor(scale_factor);
       _vars[tid].add(var_name, var);
     }
@@ -561,7 +557,7 @@ class SystemTempl : public SystemBase
     for (THREAD_ID tid = 0; tid < libMesh::n_threads(); tid++)
     {
       //FIXME: we cannot refer fetype in libMesh at this point, so we will just make a copy in MooseVariableBase.
-      MooseVariableScalar * var = new MooseVariableScalar(var_num, type, *this, _subproblem.assembly(tid), _var_kind);
+      MooseVariableScalar * var = new MooseVariableScalar(var_name, var_num, type, *this, _subproblem.assembly(tid), _var_kind);
       var->scalingFactor(scale_factor);
       _vars[tid].add(var_name, var);
     }
diff --git a/framework/include/bcs/ArrayDirichletBC.h b/framework/include/bcs/ArrayDirichletBC.h
new file mode 100644
index 000000000000..c8af06667dac
--- /dev/null
+++ b/framework/include/bcs/ArrayDirichletBC.h
@@ -0,0 +1,42 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#ifndef ARRAYDIRICHLETBC_H
+#define ARRAYDIRICHLETBC_H
+
+#include "ArrayNodalBC.h"
+
+class ArrayDirichletBC;
+
+template<>
+InputParameters validParams<ArrayDirichletBC>();
+
+/**
+ * Boundary condition of a Dirichlet type
+ *
+ * Sets the value in the node
+ */
+class ArrayDirichletBC : public ArrayNodalBC
+{
+public:
+  ArrayDirichletBC(const InputParameters & parameters);
+
+protected:
+  virtual void computeQpResidual() override;
+
+  /// The value for this BC
+  Eigen::VectorXd _value;
+};
+
+#endif /* ARRAYDIRICHLETBC_H */
diff --git a/framework/include/bcs/ArrayNodalBC.h b/framework/include/bcs/ArrayNodalBC.h
new file mode 100644
index 000000000000..7daf325ffb3c
--- /dev/null
+++ b/framework/include/bcs/ArrayNodalBC.h
@@ -0,0 +1,76 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#ifndef ARRAYNODALBC_H
+#define ARRAYNODALBC_H
+
+#include "NodalBCBase.h"
+
+// Forward declarations
+class ArrayNodalBC;
+
+// libMesh forward declarations
+namespace libMesh
+{
+template <typename T> class NumericVector;
+}
+
+template<>
+InputParameters validParams<ArrayNodalBC>();
+
+/**
+ * Base class for deriving any boundary condition that works at nodes
+ */
+class ArrayNodalBC : public NodalBCBase
+{
+public:
+  ArrayNodalBC(const InputParameters & parameters);
+
+  virtual void computeResidual(NumericVector<Number> & residual);
+  virtual void computeJacobian();
+  virtual void computeOffDiagJacobian(unsigned int jvar);
+
+protected:
+  /// The MooseVariable this BC is acting on
+  ArrayMooseVariable & _array_var;
+
+  /// current node being processed
+  const Node * & _current_node;
+
+  /// The residual values to be set in computeQpResidual()
+  Eigen::Map<Eigen::VectorXd> _residual;
+
+  /// Quadrature point index
+  unsigned int _qp;
+
+  /// Value of the unknown variable this BC is acting on
+  const ArrayVariableValue & _u;
+
+  virtual void computeQpResidual() = 0;
+
+  /**
+   * The user can override this function to compute the "on-diagonal"
+   * Jacobian contribution for this ArrayNodalBC.  If not overriden,
+   * returns 1.
+   */
+  virtual Real computeQpJacobian();
+
+  /**
+   * This is the virtual that derived classes should override for
+   * computing an off-diagonal jacobian component.
+   */
+  virtual Real computeQpOffDiagJacobian(unsigned int jvar);
+};
+
+#endif /* ARRAYNODALBC_H */
diff --git a/framework/include/bcs/NodalBC.h b/framework/include/bcs/NodalBC.h
index f136254d8c3b..bad0d55c3f95 100644
--- a/framework/include/bcs/NodalBC.h
+++ b/framework/include/bcs/NodalBC.h
@@ -15,9 +15,7 @@
 #ifndef NODALBC_H
 #define NODALBC_H
 
-#include "BoundaryCondition.h"
-#include "RandomInterface.h"
-#include "CoupleableMooseVariableDependencyIntermediateInterface.h"
+#include "NodalBCBase.h"
 
 // Forward declarations
 class NodalBC;
@@ -34,10 +32,7 @@ InputParameters validParams<NodalBC>();
 /**
  * Base class for deriving any boundary condition that works at nodes
  */
-class NodalBC :
-  public BoundaryCondition,
-  public RandomInterface,
-  public CoupleableMooseVariableDependencyIntermediateInterface
+class NodalBC : public NodalBCBase
 {
 public:
   NodalBC(const InputParameters & parameters);
diff --git a/framework/include/bcs/NodalBCBase.h b/framework/include/bcs/NodalBCBase.h
new file mode 100644
index 000000000000..7a5de469bd45
--- /dev/null
+++ b/framework/include/bcs/NodalBCBase.h
@@ -0,0 +1,50 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#ifndef NODALBCBASE_H
+#define NODALBCBASE_H
+
+#include "BoundaryCondition.h"
+#include "RandomInterface.h"
+#include "CoupleableMooseVariableDependencyIntermediateInterface.h"
+
+// Forward declarations
+class NodalBCBase;
+
+// libMesh forward declarations
+namespace libMesh
+{
+template <typename T> class NumericVector;
+}
+
+template<>
+InputParameters validParams<NodalBCBase>();
+
+/**
+ * Base class for deriving any boundary condition that works at nodes
+ */
+class NodalBCBase :
+  public BoundaryCondition,
+  public RandomInterface,
+  public CoupleableMooseVariableDependencyIntermediateInterface
+{
+public:
+  NodalBCBase(const InputParameters & parameters);
+
+  virtual void computeResidual(NumericVector<Number> & residual) = 0;
+  virtual void computeJacobian() = 0;
+  virtual void computeOffDiagJacobian(unsigned int jvar) = 0;
+};
+
+#endif /* NODALBCBASE_H */
diff --git a/framework/src/actions/AddArrayBCAction.C b/framework/src/actions/AddArrayBCAction.C
new file mode 100644
index 000000000000..66c38b3f2107
--- /dev/null
+++ b/framework/src/actions/AddArrayBCAction.C
@@ -0,0 +1,36 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "AddArrayBCAction.h"
+#include "FEProblem.h"
+#include "BoundaryCondition.h"
+
+template<>
+InputParameters validParams<AddArrayBCAction>()
+{
+  InputParameters params = validParams<MooseObjectAction>();
+  params += validParams<BoundaryCondition>();
+  return params;
+}
+
+AddArrayBCAction::AddArrayBCAction(InputParameters params) :
+    MooseObjectAction(params)
+{
+}
+
+void
+AddArrayBCAction::act()
+{
+  _problem->addBoundaryCondition(_type, _name, _moose_object_pars);
+}
diff --git a/framework/src/base/ArrayMooseVariable.C b/framework/src/base/ArrayMooseVariable.C
index 9bd45d5e3599..86be5c0d5ec3 100644
--- a/framework/src/base/ArrayMooseVariable.C
+++ b/framework/src/base/ArrayMooseVariable.C
@@ -28,8 +28,8 @@
 #include "libmesh/dense_vector.h"
 #include "libmesh/petsc_vector.h"
 
-ArrayMooseVariable::ArrayMooseVariable(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count) :
-    MooseVariableBase(0 /* FIXME: use the libmesh variable number here! */, fe_type, sys, assembly, var_kind, count),
+ArrayMooseVariable::ArrayMooseVariable(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count) :
+    MooseVariableBase(name, var_num, fe_type, sys, assembly, var_kind, count),
 
     _qrule(_assembly.qRule()),
     _qrule_face(_assembly.qRuleFace()),
@@ -38,7 +38,7 @@ ArrayMooseVariable::ArrayMooseVariable(unsigned int var_num, const FEType & fe_t
     _neighbor(_assembly.neighbor()),
 
     _mapped_values(NULL, 0),
-    _mapped_grad_phi(NULL),
+    _mapped_grad_phi(0),
 
     _need_u_old(false),
     _need_u_older(false),
@@ -205,7 +205,6 @@ void
 ArrayMooseVariable::prepareIC()
 {
   _dof_map.dof_indices(_elem, _dof_indices, _var_num);
-  _nodal_u.resize(_dof_indices.size());
 
   unsigned int nqp = _qrule->n_points();
   _u.resize(nqp);
@@ -604,7 +603,6 @@ ArrayMooseVariable::computeElemValues()
   }
 
   const_cast<PetscVector<Real> &>(current_solution).restore_array();
-
 }
 
 // void
@@ -1088,43 +1086,42 @@ ArrayMooseVariable::computeElemValues()
 //   }
 // }
 
-// void
-// ArrayMooseVariable::computeNodalValues()
-// {
-//   if (_is_defined)
-//   {
-//     const unsigned int n = _dof_indices.size();
-//     mooseAssert(n, "Defined but empty?");
-//     _nodal_u.resize(n);
-//     _sys.currentSolution()->get(_dof_indices, &_nodal_u[0]);
-//     if (_subproblem.isTransient())
-//     {
-//       _nodal_u_old.resize(n);
-//       _nodal_u_older.resize(n);
-//       _sys.solutionOld().get(_dof_indices, &_nodal_u_old[0]);
-//       _sys.solutionOlder().get(_dof_indices, &_nodal_u_older[0]);
+void
+ArrayMooseVariable::computeNodalValues()
+{
+  if (_is_defined)
+  {
+    _nodal_u.resize(1);
+    _nodal_u[0].resize(_count);
 
-//       _nodal_u_dot.resize(n);
-//       _nodal_du_dot_du.resize(n);
-//       for (unsigned int i = 0; i < n; i++)
-//       {
-//         _nodal_u_dot[i] = _sys.solutionUDot()(_dof_indices[i]);
-//         _nodal_du_dot_du[i] = _sys.duDotDu();
-//       }
-//     }
-//   }
-//   else
-//   {
-//     _nodal_u.resize(0);
-//     if (_subproblem.isTransient())
-//     {
-//       _nodal_u_old.resize(0);
-//       _nodal_u_older.resize(0);
-//       _nodal_u_dot.resize(0);
-//       _nodal_du_dot_du.resize(0);
-//     }
-//   }
-// }
+    auto & petsc_solution = dynamic_cast<const PetscVector<Number> &>(*_sys.currentSolution());
+
+    auto local_nodal_dof_index = petsc_solution.map_global_to_local_index(_nodal_dof_index);
+
+    // About the const_cast: See the explanation for the same thing in computeElemValues()
+    auto solution_values = const_cast<PetscScalar *>(petsc_solution.get_array_read());
+
+    // Ok - nodalSln() must return the same type as sln() so, for now, I'm going to do this copy here
+    // It would be much more awesome to just Eigen::Map the values directly... but to do that I would
+    // have to change the type of _u to be a Map and then I would have to compute the values into an
+    // Eigen::Vector (which I would store elsewhere) and then Map them.  That sounds like no fun...
+    // so for now I'm just doing a copy...
+    std::memcpy(_nodal_u[0].data(), solution_values + local_nodal_dof_index, _count * sizeof(Real));
+
+    const_cast<PetscVector<Real> &>(petsc_solution).restore_array();
+  }
+  else
+  {
+    _nodal_u.resize(0);
+    if (_subproblem.isTransient())
+    {
+      _nodal_u_old.resize(0);
+      _nodal_u_older.resize(0);
+      _nodal_u_dot.resize(0);
+      _nodal_du_dot_du.resize(0);
+    }
+  }
+}
 
 // void
 // ArrayMooseVariable::computeNodalNeighborValues()
diff --git a/framework/src/base/Assembly.C b/framework/src/base/Assembly.C
index aa79541ee9ac..c1a92afb2e4c 100644
--- a/framework/src/base/Assembly.C
+++ b/framework/src/base/Assembly.C
@@ -921,8 +921,6 @@ Assembly::init()
 
   auto n_array_vars = _sys.numArrayMooseVariables();
 
-  std::cout<<"n_array_vars: "<<n_array_vars<<std::endl;
-
   // two vectors: one for time residual contributions and one for non-time residual contributions
   _array_sub_Re.resize(2);
   for (auto i = decltype(_array_sub_Re.size())(0); i < _array_sub_Re.size(); i++)
@@ -1203,7 +1201,9 @@ Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
       _tmp_Re = res_block;
       _tmp_Re *= scaling_factor;
 
-      for (unsigned int i=0; i<_tmp_Re.size(); i++)
+      auto tmp_Re_size = _tmp_Re.size();
+
+      for (unsigned int i=0; i<_tmp_Re_size; i++)
       {
         cached_residual_values.push_back(_tmp_Re(i));
         cached_residual_rows.push_back(_temp_dof_indices[i]);
@@ -1211,7 +1211,9 @@ Assembly::cacheResidualBlock(std::vector<Real> & cached_residual_values,
     }
     else
     {
-      for (unsigned int i=0; i<res_block.size(); i++)
+      auto res_block_size = res_block.size();
+
+      for (unsigned int i=0; i<res_block_size; i++)
       {
         cached_residual_values.push_back(res_block(i));
         cached_residual_rows.push_back(_temp_dof_indices[i]);
diff --git a/framework/src/base/Moose.C b/framework/src/base/Moose.C
index 9fcdf91c94d9..d36e21bf1923 100644
--- a/framework/src/base/Moose.C
+++ b/framework/src/base/Moose.C
@@ -88,6 +88,9 @@
 #include "PostprocessorDirichletBC.h"
 #include "OneDEqualValueConstraintBC.h"
 
+// ArrayBCs
+#include "ArrayDirichletBC.h"
+
 // auxkernels
 #include "ConstantAux.h"
 #include "FunctionAux.h"
@@ -343,6 +346,7 @@
 
 // Actions
 #include "AddBCAction.h"
+#include "AddArrayBCAction.h"
 #include "AddDiracKernelAction.h"
 #include "AddICAction.h"
 #include "AddInitialConditionAction.h"
@@ -509,6 +513,9 @@ registerObjects(Factory & factory)
   registerBoundaryCondition(PostprocessorDirichletBC);
   registerBoundaryCondition(OneDEqualValueConstraintBC);
 
+  // ArrayBCs
+  registerArrayBoundaryCondition(ArrayDirichletBC);
+
   // dirac kernels
   registerDiracKernel(ConstantPointSource);
 
@@ -840,6 +847,7 @@ addActionTypes(Syntax & syntax)
 
   registerMooseObjectTask("add_material",                 Material,               false);
   registerMooseObjectTask("add_bc",                       BoundaryCondition,      false);
+  registerMooseObjectTask("add_array_bc",                 BoundaryCondition,      false);
   registerMooseObjectTask("add_function",                 Function,               false);
 
   registerMooseObjectTask("add_aux_kernel",               AuxKernel,              false);
@@ -988,7 +996,7 @@ addActionTypes(Syntax & syntax)
 "(add_output)"
 "(add_postprocessor)"
 "(add_vector_postprocessor)"
-"(add_aux_kernel, add_bc, add_damper, add_dirac_kernel, add_kernel, add_nodal_kernel, add_dg_kernel, add_interface_kernel, add_scalar_kernel, add_array_kernel, add_aux_scalar_kernel, add_indicator, add_marker)"
+"(add_aux_kernel, add_bc, add_array_bc, add_damper, add_dirac_kernel, add_kernel, add_nodal_kernel, add_dg_kernel, add_interface_kernel, add_scalar_kernel, add_array_kernel, add_aux_scalar_kernel, add_indicator, add_marker)"
 "(add_control)"
 "(check_output)"
 "(check_integrity)"
@@ -1078,6 +1086,7 @@ registerActions(Syntax & syntax, ActionFactory & action_factory)
   registerAction(AddDGKernelAction, "add_dg_kernel");
   registerAction(AddInterfaceKernelAction, "add_interface_kernel");
   registerAction(AddBCAction, "add_bc");
+  registerAction(AddArrayBCAction, "add_array_bc");
   registerAction(EmptyAction, "no_action");  // placeholder
   registerAction(AddPeriodicBCAction, "add_periodic_bc");
   registerAction(AddMaterialAction, "add_material");
diff --git a/framework/src/base/MooseVariable.C b/framework/src/base/MooseVariable.C
index 60653740f6f1..f7d8d8ece5a1 100644
--- a/framework/src/base/MooseVariable.C
+++ b/framework/src/base/MooseVariable.C
@@ -27,8 +27,8 @@
 #include "libmesh/quadrature.h"
 #include "libmesh/dense_vector.h"
 
-MooseVariable::MooseVariable(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind) :
-    MooseVariableBase(var_num, fe_type, sys, assembly, var_kind),
+MooseVariable::MooseVariable(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind) :
+    MooseVariableBase(name, var_num, fe_type, sys, assembly, var_kind),
 
     _qrule(_assembly.qRule()),
     _qrule_face(_assembly.qRuleFace()),
@@ -908,6 +908,7 @@ MooseVariable::computeElemValues()
     for (unsigned int qp=0; qp < nqp; qp++)
     {
       phi_local = _phi[i][qp];
+
       dphi_qp = &_grad_phi[i][qp];
 
       grad_u_qp = &_grad_u[qp];
diff --git a/framework/src/base/MooseVariableBase.C b/framework/src/base/MooseVariableBase.C
index 357e587e4dc6..86c8a1624c1b 100644
--- a/framework/src/base/MooseVariableBase.C
+++ b/framework/src/base/MooseVariableBase.C
@@ -23,7 +23,8 @@
 #include "libmesh/dof_map.h"
 
 
-MooseVariableBase::MooseVariableBase(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count) :
+MooseVariableBase::MooseVariableBase(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind, unsigned int count) :
+    _name(name),
     _var_num(var_num),
     _fe_type(fe_type),
     _var_kind(var_kind),
@@ -45,9 +46,25 @@ MooseVariableBase::~MooseVariableBase()
 const std::string &
 MooseVariableBase::name() const
 {
-  return _sys.system().variable(_var_num).name();
+  return _name;
 }
 
+const std::vector<std::string>
+MooseVariableBase::names() const
+{
+  if (count() == 1)
+    return std::vector<std::string>{_name};
+
+  // For ArrayVariables
+  std::vector<std::string> the_names(count());
+
+  for (auto i = decltype(count())(0); i < count(); i++)
+    the_names[i] = _name + "_" + std::to_string(i);
+
+  return the_names;
+}
+
+
 Order
 MooseVariableBase::order() const
 {
diff --git a/framework/src/base/MooseVariableScalar.C b/framework/src/base/MooseVariableScalar.C
index 7233c8a16307..262e562b7cd2 100644
--- a/framework/src/base/MooseVariableScalar.C
+++ b/framework/src/base/MooseVariableScalar.C
@@ -24,8 +24,8 @@
 #include "libmesh/numeric_vector.h"
 #include "libmesh/dof_map.h"
 
-MooseVariableScalar::MooseVariableScalar(unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind) :
-    MooseVariableBase(var_num, fe_type, sys, assembly, var_kind)
+MooseVariableScalar::MooseVariableScalar(const std::string & name, unsigned int var_num, const FEType & fe_type, SystemBase & sys, Assembly & assembly, Moose::VarKindType var_kind) :
+    MooseVariableBase(name, var_num, fe_type, sys, assembly, var_kind)
 {
 }
 
diff --git a/framework/src/base/NonlinearSystem.C b/framework/src/base/NonlinearSystem.C
index 78ba6d7154aa..f1c2af3fc5c9 100644
--- a/framework/src/base/NonlinearSystem.C
+++ b/framework/src/base/NonlinearSystem.C
@@ -573,7 +573,7 @@ NonlinearSystem::addBoundaryCondition(const std::string & bc_name, const std::st
   _vars[tid].addBoundaryVar(boundary_ids, &bc->variable());
 
   // Cast to the various types of BCs
-  MooseSharedPointer<NodalBC> nbc = MooseSharedNamespace::dynamic_pointer_cast<NodalBC>(bc);
+  MooseSharedPointer<NodalBCBase> nbc = MooseSharedNamespace::dynamic_pointer_cast<NodalBCBase>(bc);
   MooseSharedPointer<IntegratedBC> ibc = MooseSharedNamespace::dynamic_pointer_cast<IntegratedBC>(bc);
 
   // NodalBC
@@ -1394,7 +1394,7 @@ NonlinearSystem::computeNodalBCs(NumericVector<Number> & residual)
 
         if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
         {
-          const std::vector<MooseSharedPointer<NodalBC> > & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
+          auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
           for (const auto & nbc : bcs)
             if (nbc->shouldApply())
               nbc->computeResidual(residual);
@@ -2013,7 +2013,7 @@ NonlinearSystem::computeJacobianInternal(SparseMatrix<Number> &  jacobian)
       // safe if there are NodalBCs there to be gotten...
       if (_nodal_bcs.hasActiveBoundaryObjects(bid))
       {
-        const std::vector<MooseSharedPointer<NodalBC> > & bcs = _nodal_bcs.getActiveBoundaryObjects(bid);
+        auto & bcs = _nodal_bcs.getActiveBoundaryObjects(bid);
         for (const auto & bc : bcs)
         {
           const std::vector<MooseVariable *> & coupled_moose_vars = bc->getCoupledMooseVars();
@@ -2047,7 +2047,7 @@ NonlinearSystem::computeJacobianInternal(SparseMatrix<Number> &  jacobian)
       {
         _fe_problem.reinitNodeFace(node, boundary_id, 0);
 
-        const std::vector<MooseSharedPointer<NodalBC> > & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
+        auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
         for (const auto & bc : bcs)
         {
           // Get the set of involved MOOSE vars for this BC
@@ -2181,7 +2181,7 @@ NonlinearSystem::computeJacobianBlocks(std::vector<JacobianBlock *> & blocks)
 
         if (_nodal_bcs.hasActiveBoundaryObjects(boundary_id))
         {
-          const std::vector<MooseSharedPointer<NodalBC> > & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
+          auto & bcs = _nodal_bcs.getActiveBoundaryObjects(boundary_id);
 
           if (node->processor_id() == processor_id())
           {
diff --git a/framework/src/base/SystemBase.C b/framework/src/base/SystemBase.C
index 029cb58afee5..172cfc8ae683 100644
--- a/framework/src/base/SystemBase.C
+++ b/framework/src/base/SystemBase.C
@@ -210,12 +210,25 @@ Order
 SystemBase::getMinQuadratureOrder()
 {
   Order order = CONSTANT;
-  std::vector<MooseVariable *> vars = _vars[0].variables();
-  for (const auto & var : vars)
+
+  {
+    std::vector<MooseVariable *> vars = _vars[0].variables();
+    for (const auto & var : vars)
+    {
+      FEType fe_type = var->feType();
+      if (fe_type.default_quadrature_order() > order)
+        order = fe_type.default_quadrature_order();
+    }
+  }
+
   {
-    FEType fe_type = var->feType();
-    if (fe_type.default_quadrature_order() > order)
-      order = fe_type.default_quadrature_order();
+    std::vector<ArrayMooseVariable *> vars = _vars[0].arrayVars();
+    for (const auto & var : vars)
+    {
+      FEType fe_type = var->feType();
+      if (fe_type.default_quadrature_order() > order)
+        order = fe_type.default_quadrature_order();
+    }
   }
 
   return order;
@@ -338,13 +351,27 @@ SystemBase::reinitNode(const Node * /*node*/, THREAD_ID tid)
 void
 SystemBase::reinitNodeFace(const Node * /*node*/, BoundaryID /*bnd_id*/, THREAD_ID tid)
 {
-  const std::vector<MooseVariable *> & vars = _vars[tid].variables();
-  for (const auto & var : vars)
   {
-    if (var->isNodal())
+    auto & vars = _vars[tid].variables();
+    for (const auto & var : vars)
     {
-      var->reinitNode();
-      var->computeNodalValues();
+      if (var->isNodal())
+      {
+        var->reinitNode();
+        var->computeNodalValues();
+      }
+    }
+  }
+
+  {
+    auto & vars = _vars[tid].arrayVars();
+    for (const auto & var : vars)
+    {
+      if (var->isNodal())
+      {
+        var->reinitNode();
+        var->computeNodalValues();
+      }
     }
   }
 }
diff --git a/framework/src/bcs/ArrayDirichletBC.C b/framework/src/bcs/ArrayDirichletBC.C
new file mode 100644
index 000000000000..56a90522eadf
--- /dev/null
+++ b/framework/src/bcs/ArrayDirichletBC.C
@@ -0,0 +1,37 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "ArrayDirichletBC.h"
+
+template<>
+InputParameters validParams<ArrayDirichletBC>()
+{
+  InputParameters p = validParams<ArrayNodalBC>();
+  p.addRequiredParam<Real>("value", "Value of the BC");
+  return p;
+}
+
+
+ArrayDirichletBC::ArrayDirichletBC(const InputParameters & parameters) :
+  ArrayNodalBC(parameters),
+  _value(_array_var.count())
+{
+  _value.setConstant(getParam<Real>("value"));
+}
+
+void
+ArrayDirichletBC::computeQpResidual()
+{
+  _residual.noalias() = _u[_qp] - _value;
+}
diff --git a/framework/src/bcs/ArrayNodalBC.C b/framework/src/bcs/ArrayNodalBC.C
new file mode 100644
index 000000000000..ca13c4c0d130
--- /dev/null
+++ b/framework/src/bcs/ArrayNodalBC.C
@@ -0,0 +1,141 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "ArrayNodalBC.h"
+#include "MooseVariable.h"
+#include "Assembly.h"
+
+// libMesh
+#include "libmesh/petsc_vector.h"
+
+template<>
+InputParameters validParams<ArrayNodalBC>()
+{
+  InputParameters params = validParams<NodalBCBase>();
+
+  return params;
+}
+
+
+ArrayNodalBC::ArrayNodalBC(const InputParameters & parameters) :
+    NodalBCBase(parameters),
+    _array_var(dynamic_cast<ArrayMooseVariable &>(_var)),
+    _current_node(_array_var.node()),
+    _residual(NULL, 0),
+    _u(_array_var.nodalSln())
+{
+}
+
+void
+ArrayNodalBC::computeResidual(NumericVector<Number> & residual)
+{
+  if (_array_var.isNodalDefined())
+  {
+    auto & petsc_residual = dynamic_cast<PetscVector<Number> &>(residual);
+
+    // We know we are dealing with a locally owned node
+    // Therefore this processor also owns all of the degrees of freedom on that node
+    // Therefore we can compute and set the residual directly for all of those degrees of freedom
+    //
+    // To do that we're going to:
+    // 1. Get the global dof_index for the first variable
+    // 2. Get the local dof_index for the first variable
+    // 3. Get the raw PETSc data array
+    // 4. Wrap the part of the raw residual vector covered by these dofs in an Eigen Vector
+    // 5. Compute the residual for every variable in the array directly in place
+    // 6. ???
+    // 7. Profit!!!
+
+    // 1. This is the global dof index for the first variable
+    auto dof_idx = _array_var.nodalDofIndex();
+
+    // 2. Map it to the local dof index:
+    auto local_dof_idx = petsc_residual.map_global_to_local_index(dof_idx);
+
+    // 3. Grab the raw PETSc data array
+    PetscScalar * residual_values = const_cast<PetscScalar *>(petsc_residual.get_array());
+
+    // 4. Wrap up the the raw residual vector in an Eigen Map
+    new (&_residual) Eigen::Map<Eigen::VectorXd>(residual_values + local_dof_idx, _array_var.count());
+
+    // 5. Compute the residual directly into _residual
+    _qp = 0;
+    computeQpResidual();
+
+    // 6. ???
+    petsc_residual.restore_array();
+
+    // 7. Profit!!!
+  }
+}
+
+void
+ArrayNodalBC::computeJacobian()
+{
+  /*
+  // We call the user's computeQpJacobian() function and store the
+  // results in the _assembly object. We can't store them directly in
+  // the element stiffness matrix, as they will only be inserted after
+  // all the assembly is done.
+  if (_moose_var.isNodalDefined())
+  {
+    _qp = 0;
+    Real cached_val = computeQpJacobian();
+    dof_id_type cached_row = _moose_var.nodalDofIndex();
+
+    // Cache the user's computeQpJacobian() value for later use.
+    _fe_problem.assembly(0).cacheJacobianContribution(cached_row, cached_row, cached_val);
+
+    if (_has_diag_save_in)
+    {
+      Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
+      for (unsigned int i=0; i<_diag_save_in.size(); i++)
+        _diag_save_in[i]->sys().solution().set(_diag_save_in[i]->nodalDofIndex(), cached_val);
+    }
+  }
+  */
+}
+
+void
+ArrayNodalBC::computeOffDiagJacobian(unsigned int jvar)
+{
+  /*
+  if (jvar == _var.number())
+    computeJacobian();
+  else
+  {
+    _qp = 0;
+    Real cached_val = computeQpOffDiagJacobian(jvar);
+    dof_id_type cached_row = _moose_var.nodalDofIndex();
+    // Note: this only works for Lagrange variables...
+    dof_id_type cached_col = _current_node->dof_number(_sys.number(), jvar, 0);
+
+    // Cache the user's computeQpJacobian() value for later use.
+    _fe_problem.assembly(0).cacheJacobianContribution(cached_row, cached_col, cached_val);
+  }
+  */
+}
+
+
+Real
+ArrayNodalBC::computeQpJacobian()
+{
+  return 1.;
+}
+
+Real
+ArrayNodalBC::computeQpOffDiagJacobian(unsigned int /*jvar*/)
+{
+  return 0.;
+}
diff --git a/framework/src/bcs/NodalBC.C b/framework/src/bcs/NodalBC.C
index d060b45a1ce4..618b5300b62e 100644
--- a/framework/src/bcs/NodalBC.C
+++ b/framework/src/bcs/NodalBC.C
@@ -19,19 +19,14 @@
 template<>
 InputParameters validParams<NodalBC>()
 {
-  InputParameters params = validParams<BoundaryCondition>();
-  params += validParams<RandomInterface>();
-  params.addParam<std::vector<AuxVariableName> >("save_in", "The name of auxiliary variables to save this BC's residual contributions to.  Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)");
-  params.addParam<std::vector<AuxVariableName> >("diag_save_in", "The name of auxiliary variables to save this BC's diagonal jacobian contributions to.  Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)");
+  InputParameters params = validParams<NodalBCBase>();
 
   return params;
 }
 
 
 NodalBC::NodalBC(const InputParameters & parameters) :
-    BoundaryCondition(parameters, true), // true is for being Nodal
-    RandomInterface(parameters, _fe_problem, _tid, true),
-    CoupleableMooseVariableDependencyIntermediateInterface(this, true),
+    NodalBCBase(parameters),
     _moose_var(dynamic_cast<MooseVariable &>(*_variable)),
     _current_node(_moose_var.node()),
     _u(_moose_var.nodalSln()),
diff --git a/framework/src/bcs/NodalBCBase.C b/framework/src/bcs/NodalBCBase.C
new file mode 100644
index 000000000000..c85eb7e64903
--- /dev/null
+++ b/framework/src/bcs/NodalBCBase.C
@@ -0,0 +1,34 @@
+/****************************************************************/
+/*               DO NOT MODIFY THIS HEADER                      */
+/* MOOSE - Multiphysics Object Oriented Simulation Environment  */
+/*                                                              */
+/*           (c) 2010 Battelle Energy Alliance, LLC             */
+/*                   ALL RIGHTS RESERVED                        */
+/*                                                              */
+/*          Prepared by Battelle Energy Alliance, LLC           */
+/*            Under Contract No. DE-AC07-05ID14517              */
+/*            With the U. S. Department of Energy               */
+/*                                                              */
+/*            See COPYRIGHT for full restrictions               */
+/****************************************************************/
+
+#include "NodalBCBase.h"
+
+template<>
+InputParameters validParams<NodalBCBase>()
+{
+  InputParameters params = validParams<BoundaryCondition>();
+  params += validParams<RandomInterface>();
+
+  params.addParam<std::vector<AuxVariableName> >("save_in", "The name of auxiliary variables to save this BC's residual contributions to.  Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)");
+  params.addParam<std::vector<AuxVariableName> >("diag_save_in", "The name of auxiliary variables to save this BC's diagonal jacobian contributions to.  Everything about that variable must match everything about this variable (the type, what blocks it's on, etc.)");
+
+  return params;
+}
+
+NodalBCBase::NodalBCBase(const InputParameters & parameters) :
+    BoundaryCondition(parameters, true), // true is for being Nodal
+    RandomInterface(parameters, _fe_problem, _tid, true),
+    CoupleableMooseVariableDependencyIntermediateInterface(this, true)
+{
+}
diff --git a/framework/src/outputs/AdvancedOutput.C b/framework/src/outputs/AdvancedOutput.C
index 46b59e64e25d..b373c9234bc4 100644
--- a/framework/src/outputs/AdvancedOutput.C
+++ b/framework/src/outputs/AdvancedOutput.C
@@ -467,11 +467,17 @@ AdvancedOutput<T>::initAvailableLists()
     if (T::_problem_ptr->hasVariable(var_name))
     {
       MooseVariableBase & var = T::_problem_ptr->getVariableBase(0, var_name);
-      const FEType type = var.feType();
-      if (type.order == CONSTANT)
-        _execute_data["elemental"].available.insert(var_name);
-      else
-        _execute_data["nodal"].available.insert(var_name);
+
+      auto names = var.names();
+
+      for (auto & name : names)
+      {
+        const FEType type = var.feType();
+        if (type.order == CONSTANT)
+          _execute_data["elemental"].available.insert(name);
+        else
+          _execute_data["nodal"].available.insert(name);
+      }
     }
 
     else if (T::_problem_ptr->hasScalarVariable(var_name))
@@ -501,7 +507,6 @@ template<class T>
 void
 AdvancedOutput<T>::initShowHideLists(const std::vector<VariableName> & show, const std::vector<VariableName> & hide)
 {
-
   // Storage for user-supplied input that is unknown as a variable or postprocessor
   std::set<std::string> unknown;
 
@@ -515,12 +520,18 @@ AdvancedOutput<T>::initShowHideLists(const std::vector<VariableName> & show, con
   {
     if (T::_problem_ptr->hasVariable(var_name))
     {
-      MooseVariable & var = T::_problem_ptr->getVariable(0, var_name);
-      const FEType type = var.feType();
-      if (type.order == CONSTANT)
-        _execute_data["elemental"].show.insert(var_name);
-      else
-        _execute_data["nodal"].show.insert(var_name);
+      MooseVariableBase & var = T::_problem_ptr->getVariableBase(0, var_name);
+
+      auto names = var.names();
+
+      for (auto & name : names)
+      {
+        const FEType type = var.feType();
+        if (type.order == CONSTANT)
+          _execute_data["elemental"].show.insert(name);
+        else
+          _execute_data["nodal"].show.insert(name);
+      }
     }
     else if (T::_problem_ptr->hasScalarVariable(var_name))
       _execute_data["scalars"].show.insert(var_name);
@@ -538,11 +549,17 @@ AdvancedOutput<T>::initShowHideLists(const std::vector<VariableName> & show, con
     if (T::_problem_ptr->hasVariable(var_name))
     {
       MooseVariable & var = T::_problem_ptr->getVariable(0, var_name);
-      const FEType type = var.feType();
-      if (type.order == CONSTANT)
-        _execute_data["elemental"].hide.insert(var_name);
-      else
-        _execute_data["nodal"].hide.insert(var_name);
+
+      auto names = var.names();
+
+      for (auto & name : names)
+      {
+        const FEType type = var.feType();
+        if (type.order == CONSTANT)
+          _execute_data["elemental"].hide.insert(name);
+        else
+          _execute_data["nodal"].hide.insert(name);
+      }
     }
     else if (T::_problem_ptr->hasScalarVariable(var_name))
       _execute_data["scalars"].hide.insert(var_name);
diff --git a/framework/src/parser/MooseSyntax.C b/framework/src/parser/MooseSyntax.C
index 95ebd339a56c..7e9502cd6dcd 100644
--- a/framework/src/parser/MooseSyntax.C
+++ b/framework/src/parser/MooseSyntax.C
@@ -43,6 +43,7 @@ void associateSyntax(Syntax & syntax, ActionFactory & action_factory)
   syntax.registerActionSyntax("AddArrayKernelAction", "ArrayKernels/*", "add_array_kernel");
 
   syntax.registerActionSyntax("AddBCAction", "BCs/*", "add_bc");
+  syntax.registerActionSyntax("AddArrayBCAction", "ArrayBCs/*", "add_array_bc");
 
   syntax.registerActionSyntax("CreateProblemAction", "Problem");
   syntax.registerActionSyntax("DynamicObjectRegistrationAction", "Problem");
diff --git a/test/tests/arraykernels/simple_diffusion/gold/simple_diffusion_out.e b/test/tests/arraykernels/simple_diffusion/gold/simple_diffusion_out.e
new file mode 100644
index 000000000000..5de69d00460f
Binary files /dev/null and b/test/tests/arraykernels/simple_diffusion/gold/simple_diffusion_out.e differ
diff --git a/test/tests/arraykernels/simple_diffusion/simple_diffusion.i b/test/tests/arraykernels/simple_diffusion/simple_diffusion.i
index 810df9edd7ef..1ea8ccb66cd0 100644
--- a/test/tests/arraykernels/simple_diffusion/simple_diffusion.i
+++ b/test/tests/arraykernels/simple_diffusion/simple_diffusion.i
@@ -1,13 +1,14 @@
 [Mesh]
   type = GeneratedMesh
   dim = 2
-  nx = 10
-  ny = 10
+  nx = 20
+  ny = 20
+  nz = 20
 []
 
 [Variables]
   [./u]
-    count = 3
+    count = 2
   [../]
 []
 
@@ -18,16 +19,15 @@
   [../]
 []
 
-[BCs]
-  active = ''
+[ArrayBCs]
   [./left]
-    type = DirichletBC
+    type = ArrayDirichletBC
     variable = u
     boundary = left
     value = 0
   [../]
   [./right]
-    type = DirichletBC
+    type = ArrayDirichletBC
     variable = u
     boundary = right
     value = 1
@@ -44,12 +44,12 @@
 
   # Preconditioned JFNK (default)
   solve_type = 'JFNK'
-
-
-  petsc_options_iname = '-pc_type -pc_hypre_type'
-  petsc_options_value = 'hypre boomeramg'
 []
 
 [Outputs]
   exodus = true
+  [./console]
+    type = Console
+    perf_log = true
+  [../]
 []
diff --git a/test/tests/arraykernels/simple_diffusion/standard.i b/test/tests/arraykernels/simple_diffusion/standard.i
new file mode 100644
index 000000000000..8396c5deba3c
--- /dev/null
+++ b/test/tests/arraykernels/simple_diffusion/standard.i
@@ -0,0 +1,28 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 20
+  ny = 20
+  nz = 20
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = JFNK
+[]
+
+[Outputs]
+  exodus = true
+  [./console]
+    type = Console
+    perf_log = true
+  [../]
+[]
+
+[Testing]
+  [./LotsOfDiffusion]
+    [./vars]
+      number = 2
+    [../]
+  [../]
+[]